{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "离散概率分布的表示函数称为概率质量函数\n",
    "连续概率分布的表示函数称为概率密度函数\n",
    "概率质量函数的输出是概率，概率密度函数曲线下面积表示概率"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "预估允许的最大标准差c = [31.21216584]\n"
     ]
    },
    {
     "data": {
      "text/plain": "3.1213862058843134e-05"
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#4.1\n",
    "\n",
    "#X~N(160，x**2)\n",
    "from scipy.stats import norm\n",
    "from scipy.optimize import fsolve\n",
    "f = lambda c: norm.cdf(200,160,c)-norm.cdf(120,160,c)-0.8\n",
    "print('预估允许的最大标准差c =',fsolve(f,8)) #如果给与初值0或其它不合适的初值，则会报runtime error(说明赋予的初值不合适)，并且也返回原值，所以在返回值和赋予初值相同时，需要留意！！！！\n",
    "# f(31.21)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.0, 1.2, 2.4, 3.6, 4.8, 6.0, 7.2, 8.4, 9.6, 10.8, 12.0, 13.2, 14.4, 15.6, 16.8, 18.0, 19.2, 20.4, 21.6, 22.8, 24.0, 25.2, 26.4, 27.6, 28.8, 30.0, 31.2, 32.4, 33.6, 34.8, 36.0, 37.2, 38.4, 39.6, 40.8, 42.0, 43.2, 44.4, 45.6, 46.8, 48.0, 49.2, 50.4, 51.6, 52.8, 54.0, 55.2, 56.4, 57.6, 58.8, 60.0, 61.2, 62.4, 63.6, 64.8, 66.0, 67.2, 68.4, 69.6, 70.8, 72.0, 73.2, 74.4, 75.6, 76.8, 78.0, 79.2, 80.4, 81.6, 82.8, 84.0, 85.2, 86.4, 87.6, 88.8, 90.0, 91.2, 92.4, 93.6, 94.8, 96.0, 97.2, 98.4, 99.6, 100.8, 102.0, 103.2, 104.4, 105.6, 106.8, 108.0, 109.2, 110.4, 111.6, 112.8, 114.0, 115.2, 116.4, 117.6, 118.79, 119.99, 121.19, 122.39, 123.59, 124.79, 125.99, 127.19, 128.39, 129.58, 130.78, 131.98, 133.18, 134.38, 135.57, 136.77, 137.97, 139.16, 140.36, 141.55, 142.75, 143.94, 145.13, 146.33, 147.52, 148.71, 149.9, 151.09, 152.27, 153.46, 154.65, 155.83, 157.01, 158.19, 159.37, 160.55, 161.73, 162.9, 164.07, 165.24, 166.4, 167.57, 168.73, 169.88, 171.04, 172.19, 173.33, 174.47, 175.61, 176.74, 177.87, 178.99, 180.11, 181.22, 182.32, 183.42, 184.51, 185.59, 186.66, 187.73, 188.79, 189.84, 190.88, 191.91, 192.93, 193.94, 194.94, 195.93, 196.9, 197.86, 198.81, 199.75, 200.68, 201.58, 202.48, 203.36, 204.22, 205.07, 205.9, 206.72, 207.51, 208.29, 209.06, 209.8, 210.52, 211.23, 211.91, 212.58, 213.22, 213.85, 214.45, 215.04, 215.6, 216.14, 216.65, 217.15, 217.62, 218.08, 218.5, 218.91, 219.29, 219.66, 219.99, 220.31, 220.6, 220.88, 221.12, 221.35, 221.55, 221.74, 221.9, 222.04, 222.15, 222.25, 222.32, 222.38, 222.41, 222.43, 222.42, 222.4, 222.36, 222.29, 222.21, 222.12, 222.0, 221.87, 221.72, 221.56, 221.38, 221.18, 220.98, 220.75, 220.51, 220.26, 220.0, 219.73, 219.44, 219.14, 218.83, 218.51, 218.18, 217.84, 217.49, 217.13, 216.76, 216.39, 216.01, 215.62, 215.22, 214.82, 214.41, 213.99, 213.57, 213.14, 212.71, 212.27, 211.83, 211.39, 210.94, 210.48, 210.03, 209.57, 209.1, 208.64, 208.17, 207.7, 207.23, 206.75, 206.27, 205.79, 205.31, 204.83, 204.35, 203.86, 203.37, 202.89, 202.4, 201.91, 201.42, 200.93, 200.43, 199.94, 199.45, 198.95, 198.46, 197.96, 197.47, 196.97, 196.47, 195.98, 195.48, 194.98, 194.48, 193.98, 193.49, 192.99, 192.49, 191.99, 191.49, 190.99, 190.49, 189.99, 189.49, 189.0, 188.5, 188.0, 187.5, 187.0, 186.5, 186.0, 185.5, 185.0, 184.5, 184.0, 183.5, 183.0, 182.5, 182.0, 181.5, 181.0, 180.5, 180.0, 179.5, 179.0, 178.5, 178.0, 177.5, 177.0, 176.5, 176.0, 175.5, 175.0, 174.5, 174.0, 173.5, 173.0, 172.5, 172.0, 171.5, 171.0, 170.5, 170.0, 169.5, 169.0, 168.5, 168.0, 167.5, 167.0, 166.5, 166.0, 165.5, 165.0, 164.5, 164.0, 163.5, 163.0, 162.5, 162.0, 161.5, 161.0, 160.5, 160.0, 159.5, 159.0, 158.5, 158.0, 157.5, 157.0, 156.5, 156.0, 155.5, 155.0, 154.5, 154.0, 153.5, 153.0, 152.5, 152.0, 151.5, 151.0, 150.5, 150.0, 149.5, 149.0, 148.5, 148.0, 147.5, 147.0, 146.5, 146.0, 145.5, 145.0, 144.5, 144.0, 143.5, 143.0, 142.5, 142.0, 141.5, 141.0, 140.5]\n",
      "222.43\n",
      "31.920000000000005\n",
      "216\n",
      "216 0.8674284731834105\n",
      "(21.6,22.24)\n",
      "21.6\n",
      "21 22.203574034386\n",
      "22 22.229440654350366\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAETCAYAAAA/NdFSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABn60lEQVR4nO2dd3hURReH30nvCYFQQui9996bFOmKig35VJBmRQFBBemIFQFFsIEFEKVLUVoA6UWK1NCSEEhvkLY73x93AyEkIQm7ubvJvM+TJ7v3zp357c1mzp0zZ84IKSUKhUKhKHrY6S1AoVAoFPqgDIBCoVAUUZQBUCgUiiKKMgAKhUJRRFEGQKFQKIooygAoFGZECGEnhBA6te0jhHAVQjjq0b7C9lAGoAgihCgphOjxENeXF0K8k8P554QQc/NQX00hhFsWx4UQoknGDlUIMVoI4SKE2CaEqCeEeFsI4SWE+EoI0T7vnyZLPa2FEJNyWdbHpHGYEGIlEAG8kOF8cyHEngzv6wohPE2vewohnjeHZhNjgbeB14UQXwohsv3/FkK8J4T4MDeVCiFqmn5XFkJ0zkV5eyHETiFEhdwKV+iDg94CFLpQFvhZCDEUcAQmA0mmc6WBG4ATsBR4BzhpOlcZ6AJ0M53PjhTTD0KIKcBQICpD22OklL9mKP8L8CJwRAjxOtACKAMEmLS8B2wzlXUAJgFpgDvwlJTyIyFEF2BOVmKEEFsBL+B2plOuQJyUslum4zVN7WdVlx3gKKVMNh16H6hjqusy0EhKeSXTvUjN8L4N8KMQohtwGNgrhPgd2AVI048bMFJKuSMrDTlw23T9x8CPQFtTvVmRbCqbI0KIXsBHQoi6pvKLhBD1pZS3crisB+CW6T5krvcI2t8yJYd6ykkpSz1IoyL/KANQBJFSHhVCPAmMB7pKKX9LPyeECJZSNsnw/mXgCcAZ+BQwAM8DvkKIx03F7IGOQDQgACMghRDOaB31+1LK7031zQfiMtRfD/gXOGUaBfQEPgA8gWeklC9kKOsChKJ1RMVMbe4yPfl7SCmDTOWcM3TQkHMnk9W5ssCjQogWWZyzAy4BjwFIKd80tfk64JPe6QkhngL6AB9lvFhK+bUQ4ijgAjwFNJRSJgCNM3zO73nA/6bJ4H2G1pGnUxLtbzEQzbB/B1QxlXcGUuTdlZ/GTPUJwCnjfTMZu8nAu6brLgkhNpjaHZahXCs0I56A9veuDlwXQhzL2ITpM40yGbZUYKCU8rIQojWake+Vrk8I4YBmUBUWRBmAIoqUcqsQ4i/54KXgRqAX0ND0vh1wVUrZMr2AEOK0qVw74HO0ztsNrQNcZypTDBgEdAYWZ6j/fbR//qFoo4tUNDeKRxZanID2QHOgKfAf2shiBGAnhDgElAfihRANpZTxpusOmepNylSfC1pHmZlKwOtSytVZ35I7n9sNrZNNRRs1OAkhqpl0ruTeJ3+EEB2AllLK2UKI0miuoptCiLfQOk4BZBwZ5YQ7cDiTgewLdJNSjhFC2EspDRnK/w24CSGMgDdQ1XTN82gG2Q7NeDTIcM1YIFZKuSbDsYnAPiHEx8BYqfEPUNFU32Dgf1LKbkKIxcCr2YwWUk3lvYFv0EYvB03uponAfNM9UVgQZQCKGEKIZ9FcKALNTTD9AZdItH/WWDSXTEVglhBinZSyj6mMA2CQUu4EGppGBg2llJOEEONNZRLROp0PpJRHTVpaoz3xuwK1gEfQDEhGvRmfTG8DQWhP6ElADNqoozXwpZRyqhDiV+BjKWW8EKI/8JZJe+bOPx1HIUQgWod/2HSsKrl7+kwGppnuz1fARWAu2mipUhblzwELhRCNgDFoo6lf0Dr9E0BttHu5KBdtZ2W4Q9FGAQDfCSFWSinXAUgp26YXEkJ8iXYvBbBZSjk2c0WmEcbraC6ru41KmWCaB9iKNvp6SUp51nRNRWAqmosQoFM2OtOxR3tAiAICgX3AS8A81PxkgaBuchFDSrlMSlkTzTWRkx8/HWEql+62+QEoheaySXcVOXD/01olIcQstO/YbOA6Wgf/qcnFAXAQqAucB2ZKKS+ajv+K1pH2AvYAG0zHA9BGB+8AIWiuoqtoT/HpWsqjdcRIKVdLKdsBYaZjZ9AMRrTp9UXgppSyXYbOHzRXxmIhxKEsfkKFEJNN9RuklCdMn78tmkunopTyGFl0fFLK60AHYG+G85K77hgHMrlmckAA/YUQZ4QQ14QQy033pKppkrkvcPq+i7ROuhPwM/AH0FsI4ZepjCeaEXoP+EcIESuEiBNCBAshgtGM4wLTPTSYrimD9neqAPxhcv8EoI0WjgkhooQQvTPJMaAZwqeAwWh/90+llLm9B4qHRI0Aija5yQRoB/ijdeCgPbW9CLwGLBFCjEH7HqUKIbqbjjdF6xyWormOpgAvSimbmDr/tQBSylTTZGgzYLsQ4jlTG0+hGZSvpJQZo5V8gTdM9dmhPT2moo0AdgghfABPKWUU99IYbQSQgjZ6AK2zdCILIyil7JXdzRBCLEB7kk9/b4/29L8FzaB8KoRIROug70NKGS6EmAcUz6DtSzTDtRfNuOUGN+BXKeUrppFOVynldVNH/BbwSwaDmlHrt2hzPzXQ7sdUYIUQopuUMs2kMV4IUUtKmYL2N54DhEgpPzfVswn4T0r5tel9A+B3tLmBcVLKhqbjF4AWUsok098943wFQDm0OQbQDNIeYLgQYhEwOpf3QfEQqBFAEUcI4SCE8M2hiANaZ3EOOIXmI+4KLENzc/QD7KWUqcAFYBYwCtglpdxkquMWcEAIMRpoBaw3tS2A4WiRPvWAjFEj18jkRjE9pbdFc/1UR4v6OSulvI3WAa3hbrRQ5s+wGu0J87Dp51fTsTv/A0KIckKIMCHEOdOTdcafy6aJ3Qrc6x76GM2YbAEi0eY5TqEZlvsesIQWo78fbf4BKeVBKWULtHmVplLKrPRnRRnujsq80EZGmOoewd2ONb1dO7T5itPpbiFT+z+hjRQ2ZvwemDr/dDqY6k2nLBCc4X0s2kTxvAdozvzAcQ0toKArcBaYgRa11ATY8YC6FGZAjQCKNp5onfEqtIm4rGiB5p89AqxAe7LvDBwDqpl+ADA9cV7MEB2UkSloHee09M7FNAHd1NQ5NeDuBKQH2qTqRSFEM+BMhgldf7QOfJOpzPum47+jGZ+pGRs1GZm30PzyoLl/4G4Hlj5HgZTyGloYbLYILaz1sul1Q6A70BJtEhuTSwhTtMzOLOp7Gu3eJZnKfYTJzy609QIncmo/A02AQNMkakO0e+WFZjRPSylvCCE6oRluOzSDHQzECCFOAD5oE+dPoBnNI8AZIUQLKeWlDJ833cDvy9B2GTSjh+kzX+auUbzvodIUgVSce92EwnRtiBBiHFpU1Spgp5QyzhQFpMuCuqKEMgBFF180/+s7UspvhBABaJ2re6Zys9EiMs6iuVxaSymPmM6dFULc5P5/VGfu9WsD/A/t6e4lIcRZKeVKoS0UWoD2VH0YzTVU3PR7Mpq74h3ghhAiUEq5HG2UsRzNoEwE6gshBqJNbI9Fi1MfLKVMf2KtjDZRm74GIN0FVMv021EIsV9KGZvdjTJF+3RB6/hKYHKHSSmPCSEamVwc99wD02SpE9rEbnpoowuawRqYXrWU8u0s2vslwzVhQIVM4ZmuaPMpH6CNZIqjGcTtwNdocwNvohnvv9FccT8A35mMbvrkvLOUckqGepdk0fkvAPqb3nuhzdkkZQqzzYhzhtcOaPMz/6KNVv7NcO5O9JUpKuo74Bm00eapDNcrLIi6wUUXe+BJKeUfpvdNgS/IEIFi8s97oEXYGE1PqxXRJlUx+f8no/l+0695F3gZ0xMx2j96H7Snvz5onf1yIYQBzQUzE9gjpZSmDnMe8KaU8ldTp/oa2mTwDFNHvAHtCXqglDJMCPEMMBJtQdg5IcRhYKrJL+4BfI/mHkl3PyRk+i2A9UKIx6SUN7O6UVLKW6Yn5SRTuzLDufToIkfun0/oabo/6auiawI7pJTHhRBlubezTL9/v6EZmvRO8IcsOttBpnt2EegptPDSH9BCLteaRh9b0dYAjDe5yL7NVMd9eqWU5zPo+AzNNdM3wwT5K2gum3cz686Af6Y2HIDqJhdhRhzQ3E73rcMQQryfoYzCgogHh4ErijJCCJHdWgHTpKIx43khhGPmf/bMdQgh7LKL9BBCeEkpMy4UcwX8M09oFhWEEG2llLuzOO6U7kozudCKSynDM5x3BEpJKYMzX5vLdj3RnvQzd9xmQQhRHIiR965VUBQwygAoFApFEUVFASkUCkURRRkAhUKhKKLY1CRLiRIlZMWKFfWWoVAoFDbF4cOHI6SUfpmP25QBqFixIocOHdJbhkKhUNgUQogsU3MrF5BCYSVcvnxZbwlmISgoSG8JilyiDIBCYQXMnj2bI0eOPLigDbBw4UJ27cpuHxqFNaEMgEKhM5cvX+batWt06dKFnj170q1bNwYMGEBKirZG6saNG7Rr1y5XdfXt25ejR49me/7q1at07NiRzp07M2zYMKSUWR7L7bXpnDx5kkceeQSAmTNn8vnnn2M0qqSe1o5NzQFkRWpqKsHBwSQlZZfuXZGOi4sLAQEBODqqPcOtiaVLlzJ69Gh++ukn3nzzTbp168aIESPYtGkT7dq1Y8iQISQmJj6wnp9++onKlSvTqFGjbMt8/fXXLFy4kFq1atGzZ09OnDjB8uXL7ztWv379XF1bv359pJS8+eabdwyWg4MDXbt2Zc+ePbk2XAp9sHkDEBwcjKenJxUrViRTOhZFBqSUREZGEhwcTKVKWe1VotCLixcvUrNmTWrWrHnnWHh4OCVLlsTe3p7ly5fTr1+/HOuIiorirbfeYsSIEWzfvp1OnTplWW769Lv7/0RGRlKiRIksj+X2WoDvvvuOTp06sXnz5jvnW7ZsSWBgoDIAVo7NG4CkpCTV+ecCIQTFixcnPDz8wYUVFmX10RA+2nyW0Jjb+Pu44hx1746J//zzD9HR0bRs2TKbGu7n008/ZdCgQQwfPpwJEyYQHx9P3759sy2/fPly6tSpg7+/f47HHnRtZGQky5YtY/PmzfcYAFdXV27fvp1DLQprwKJzAEKIUkLbbi+7845CiPVCiL1CiP89RDv5vbRIoe6T/qw+GsKE308QEnMbCYTE3OZoSCK/7NH2mImKimLMmDF8+23m3G05c/ToUUaNGkXp0qV54okn2LFjR7Zlg4KCmDt3Lp999lmOx3Jz7fjx45k5c+Z9bsVLly5Rrly5PH0GRQ7Eh1mkWouNAIS2CfgP3J9eOCNjgENSyslCiN+FtodpfA7lCwXBwcH4+/tjZ2ce+7tnzx4aN26Mq6urWepT5A4pJSsPBfNNYBCXIhKpVMKdl9tXZlCTACSSq3FXORV5itORp7kcd5nwW+GcCQ/BvvItPJBaHlJphzExiXeWPM7O+NZse38bg0cPJt49nlRDKo72uZuvqVq1KkFBQdSsWZNDhw5RoUKFLMtFR0czePBgvv32W7y9vbM9lttrd+7cyfnzWhLRY8eOMWnSJKZNm8bq1avvcRkp8smtKNg0Ac79CaMOgGeO21XkGYslgzPlDhfAGillx2zKrEVLV3taCDEWOCyl3J6pzDBgGED58uWbXLly73qG//77j1q1aqEnaWlpODhkbUuzOle3bl22bdtGyZIls7wG4Pvvv+fxxx9n9+7dCCHo3r17luUSEhKoV68eJ06cwMPD44FareF+FQaMRsn43/9lxaFgGpTzoWVlX/ZevM7pmMPUqHSVWw4niEyKBMDZ3plK3pXwc/Vj26nbSIMbUgpAIIQB7BK4/uM2ytb05Pya8ziX07JEl+5amg59OtCkVBO+H/U9+3fvx07YsW3bNk6fPs3o0Xd3TQwNDeWll14iNjYWNzc3fv/9d86dO8fWrVsZP/7OnjeMGzeOH3/8kRo1agAwZcoUNm7ceN8xg8FwXxtZXduhQ4c75zt27MiOHTvYvXs3mzZtYtq0aZa5+UUBKeHUH7DxbUiKgXZvaT8O92UQzxVCiMNSyqb3Hbd0NlAhxI4cDMDfaPnVY00dfZyU8tfs6mratKnMvBLYGjq04cOHc+bMGYQQxMTEEBUVReXKlQGoUKECw4cPZ8yYMXh6emI0Gjl69Cj169fH3t4eOzs7YmJi2LBhA2XLanuVXLlyhYEDB3Lo0CFu3rxJ7969CQwMxMXFBYCNGzcyc+ZM7O3tiY2NJSQkhNq1a9/RYzQaee2113jsscfu02oN96swMHfzWb7cfoHRnarSv7kDv51fydoLa4lPjUcanKnk3oShjbtTp3gdKvtUxtFOe5JvM2sbITH3+8ZLiATerGfgyaeeJCQhhP+i/uPYzWMcuXmEM1FnMEojfq5+dCjXgS7lu9CyTEsc7HIewEdFRbFmzRqGDh2aYzlz88knnzBmzBgVbZZf4q7Dhrfg7AbwbwR9v4TSdR+qSms1AGuA4aaNPd4EwqSUP2dXl7UagIzs2LGDTZs2MWvWrCzPv/fee1SqVImqVauyZcuWLJ+S+vTpw+jRo+889X/++efs2rWLX3/9FUdHR4xGI0ajkdTUVFq1asWOHTv4+++/s+zwM2Nt98sWOXo1mgEL9tK1USKi2Fb2Xd+Hg50D3Sp0o1+Vfqza48IfR8NYObwVTSveu91y+hzA7dS7afBdHe2ZObAe/RuVzdwUADFJMQSGBLL92nb2hOzhVtotSriWoHfl3vSt0pdqxapleV1ycjL29vbZjk4VVoaUcHQpbJ4EhmToNBFajgT7h//7ZWcA9P5mpG/y/RvafrD7ci6eM1PWneJ0aNyDC+aB2v5efNCnzgPLvfrqq3zxxRf3HTcYDNjb2wNa9MTmzZt5+eWXad++PYsXL2bVqlX3dNyffvopxYoVu7OoBuC1114jJCSENm3asHjxYurXr4+dnR1vv/02w4YNw8fHh3nz5uXKACgeDqNRMmH9FrwrrWB/0jl8o315rfFrDKg6gOKuxQFo2C+N/UGxTPzjJBtfa4e93d3J9/ROPmMU0Nvda2Tb+QP4uPjQp0of+lTpQ4ohhcDgQNZcXMOy08v4/tT31CtRj6drPU33Ct3vmTNwds6fu0ChA1GXYN2rcGkXVGgLfb+A4lUs3myBGQAhRGegtpTyywyHf0DbFq4d2t6p+7O82AbIGHXx888/s2/fPqSUdO3alZdeeonJkycTHBzMli1bePXVV6lduzaLFy9m+PDhfP3118ydOxdvb2+2bdtG9erVqVu3LhEREVSqVAmDQXtaHD169J0FQd988w0LFiygRo0arFixguPHj9OxY0cSEhJo0aIF8+fP1+M2FGpik2N5Y/Nsgt3W4+bgwahGY3mixhO4Otw7+e7u7MC4njV59ZejbDhxnb4N7g2r7N+obI4dfk442TvRpUIXulToQuTtSDZe2siKsyuYEDiBjw99zBM1nmBQ9UGUcM06ll9hZRgNsP9r2DYVhD30/hQavwBmChB5IFJKXX/Q9hB9AvB+UNkmTZrIzJw+ffq+Y3rQoEEDKaWU27dvl+PGjbvnXHR0tFyyZIk0Go1SSilDQ0Pl3r17pZRSGgwGuX//fpmWlnbPNatWrZJTpkyRUkp5+fJlOXDgwDvnli1bJlu3bi0nTZok//zzTymllF26dJFSSnn06FH52muvZavTWu6XrbH18lbZ7pd2su539WTrRSNl9O3oHMsbDEb5yCc7Zae522VqmsGi2gxGgwwMDpSvbH1F1v2+rmz0YyM5ee9keTXuqkXbVTwkN05LuaizlB94SblskJQxwRZrCi3a8r4+VfdcQFLKUCnlCillrN5aLIWPjw+7du3i5s2btGzZEmdnZyZNmgTAunXr+Oqrr+64idJZtmzZnTmA0NDQe2KqBwwYwNatW3MM2VOYh/iUeCbunsgbO97A3b4EiZfGMLHlJHxcfHK8zs5O8Ea3agSFJ7LxpGViuO+0JexoW7YtC7suZG3/tQyoOoA1F9bQ548+TNw9kaBYlZ3TqkhLgZ1z4Kt2EBUEAxfD08vBO3+jwodBdwNg6xgMhmyTZwGkpKQQHh7OxYsXKVWqFAkJCfj6+rJp0ybOnDlDcHAwERER92RP/OKLL3BwcKBFixYAnDt37h4D4ObmhpubG8nJyffoOH36NG+//Tbly5e3wCctepwIP8Hjax9nfdB6htcfjmvEG5R1r0L3OrmLxX6kdmnK+7rx074sU7FbhErelXiv1XtsemwTT9d6mi2Xt9B/dX/G7hzLuehzBaZDkQ0hh2FRR9g+HWr3g9EHof4g0GmRpt6TwDbP77//zhdffIGPjw8dO3a8czz9dXJyMl9++SVvvPEGoC3Yadu2LQ4ODly/fp0PPviAX3/9FTc3N6SUDBw4kOLFi7Ns2TJAi7XeunUrP/74431tP/3003h5eQEQFxdHlSpVmDNnDg0bNrToZy4K/HbuN2bsn4Gfqx8/9vwR+5SKzL2ym/d7175nUjcn7OwEg5uXZ/amM1y4GU/Vkp4WVn2Xkm4leafZO7xY90WW/beMX878wpbLW3i08qOMbDiScp5qlW6BknILdsyAf+aDR2kY/CvU6Km3KsuHgZoTWwgDfVgSExNxd89p8fTDUdjul7lJMaQwY/8MVp1fRWv/1sxuNxsfFx8+XHeapfsuc3BiV3zcnHJdX0RCMq1m/s1zLSvyfp/aD77AQsQmx7Lk5BJ+/u9nDNLA49UeZ3iD4WqyuCC4FKhF+EQFQZMXoNuH4FKw7tvswkCVC8jKsGTnr8iZuJQ4Rvw1glXnV/FyvZdZ0GUBPi4+pBqMrD0eQpeapfLU+QOU8HDmkTql+f1oMKkG/fLjezt782aTN9kwYAMDqg5g5bmV9Pq9F58f+Zy4FPOGTitMJMXCutfhh95ajP+QddDn8wLv/HNCGQCFAghLDGPIn0M4cvMIM9rO4NXGr2Jvp03MB54PJyIhhceaBOSr7n4N/Im5lcrei5HmlJwvSrmX4v1W77O2/1o6luvI4hOL6bmqJ0tOLCEpTe2pYTbObYb5LeHID9BqNIzYC5Xa663qPpQBUBR5zkef59mNz3I98ToLuy6kT5U+95z/42goxdwc6VDdL1/1t6/uh6ezA+uPh5pDrlko71WeOe3nsLLPShr4NeCzI5/RZ3Uf1getxyjVTl75JjECVr0EPz8Brj7w4l/QfTo4uemtLEuUAVAUaU6En2DIn0OQUvJDjx9oWebeHPzJaQa2n7lJ9zqlcXLI37+Li6M93WqXYvOpMFLSrKtzrelbkwVdF/Bt928p5lyMCYETeGbDMxy+cVhvabaFlHDiN5jfHE6tho7vwrCdENBEb2U5ogyAoshy7OYxhm0dhrezN0t7LaWGb437yvxzMZKE5DQeqVPqodrq3aAMcUlp7L5gnRvyNCvdjF97/8qMtjO4efsmL2x6gTe2v8HVuKt6S7N+YkPgl6dg1YtQrCIM3wUdx4FD3uaL9EAZAAvw2muvceTIEYYMGUJMTMx95+fNm8fixYsLXpjiDsduHuOVv17B18WX73p8h79H1rtgbTl9A3cne1pXebhomTZVS+DuZM/W0zcfqh5LYifs6FOlD+sHrGdUw1HsCd1DvzX9mHNwDrHJhXadZv4xGuHQd7CgJQTthO4z4MWtUEq/aK+8ogyABUjfoP6FF17ghx9+QEpJWlranfOOjo73ZGjMeE5heY7cOMLwrcPxc/Xj2+7fUto964VdRqNk6+kbdKxREhdH+yzL5BZnB3vaVivBjrM3c1w4aA24OrjySoNX2DBgA32r9GXZ6WX0+r0XS08vJdWQqrc86yDyIvzYF9a/DmUawMi90GoU2D3c96SgKVwLwf4cD2EnzFtn6XrQM+vUzukMHTqUCxcu3AnhvHjxIseOHcPHx4fk5GQGDRrE4MGDcXR0JCkpiT179gDapi8ODg6kpqayfPlySpc2724/ivs5FXGKEX+NoKRbSZZ0X0JJt+w35TkREkt4fDLdaj+c+yedLjVLsfnUDf67Hk9tfy+z1GlJ/Nz8mNJ6Ck/XfJq5h+Yy5+AcVpxdwdvN3qZ9gPVFtBQIRgPsWwDbpoO9I/T5Aho/r9tK3oelcBkAnRBC8M0331CzZk0A5s+fT0BAAHXr1uX48eP4+/uzc+fOO+eklBiNRoYPH87gwYP1lF6kCIoNYsRfIyjmUozFjyzOsfMHLfwToG018yyW6lhTiyLafvamTRiAdGr41mBRt0UEhgTy0cGPGPX3KNqWbcvbzd6msndlveUVHDdOwZrREHoEavSCRz8Gr6xdh7ZC4TIAD3hStxRCCJ555hmklDRq1Ig+ffoQFBTEgQMHqFq16p1ykZGRrFix4k6nv3DhQnr27ImPj48uuosSYYlhDN86HCEEX3f7mlLuD36qDzwfQe0yXpTwME9e/ZKeLtQr6822MzcZ1anqgy+wIoQQtA9oT6syrfj5zM98dfwrHlvzGINrDeaVBq/g5WQ7Bi3PpCVD4Mfaj4sPPP4t1Blos0/9GVFzAGbg9u3b/P7773z99df4+fnRoEED/vzzT/bs2cPzzz8PaDmBnn32WaZOnYqDgwMeHh6MHTuWAQMGEBERofMnKNxEJ0UzbOswElIS+KrrV1TwynrD9IwkJqdx5Go07cz09J9Oxxp+HL0aTVySbfrSHe0dGVJnCOsHrKdf1X4sO72M3r/3ZuW5lRiMhgdXYGsEH4KvO8DO2VD3MW1j9rqPFYrOH5QBMAtXr16lePHixMTE4OvrS8WKFTl+/DjvvvsuAEFBQXTp0oWmTZsyceJE5s2bx8cff8ycOXPo0KEDzZo14+jRozp/isJJYmoiI/4aQWhCKPM6z6NW8dzlQTpwKYpUgzSb+yed1lVKYJRw8FKUWestaIq7Fmdy68ks772cSt6V+PCfD3ly/ZMcDDuotzTzkJIIm96FxV0hOQ6eXgEDF4F7cb2VmZXC5QLSgbi4OOLi4ti/fz9vvfUWU6dOvbPl47p167h58yZxcXG89tprDBo0CNDmAVxcXHjxxRcBGDx4MNWqZb2vqyL/pBnTeGvnW5yJOsNnnT6jaen7cmFlS+D5CJwc7GiWaU/fh6VReR+cHOzYezGSLrXMM7msJ7WK1+L7Ht+z+cpmPjn0Cf/b/D8eqfAIbzV9K9vQWqsnaKeWvC36MjR9EbpOBpfC6eJSI4CHZOHChTz//PO0aNGCb775hrlz59KpUydWrlxJVFQUr732GpUrV77T+YOWuz8lJeXO+xo1amBXUFvAFRGklMw6MIs9IXuY2HIiHct1zNP1uy+E07yi70OHf2bGxdGephWKWUVeIHMhhKBHxR6s6b+GkQ1Hsit4F31X9+XLo19yK/WW3vJyz+0YWDtGC+8U9vDCBuj9SaHt/AH9t4TMy481bgmZkpIiU1JSpJRSGo1GGRoaes/5s2fP6iErW/S+XwXF0lNLZd3v68q5B+fm+dqw2Nuywrj1cuGOCxZQJuUXf52TFcatl1EJyRapX2+uJ1yX7+x8R9b9vq7svKKzXH9x/Z3tUK2W/zZIObeGlJN9pNzynpQpt/RWZFaw1i0hbR1HR0ccHR0B7UmoTJky95yvXr26HrKKNDuu7WDOwTl0Kd+FN5q8kefr91zQJuXbVrVMrvzWVTU/8r6gwjMKyEhp99LMbj+bH3v+SAnXEowPHM/zfz7PqYhTeku7n4RwWDkUfh0MbsXhpb+1fP2OrnorKxCUAVAUKv6L/I93dr1D7eK1mdluJnYi71/x/UFReLs6UruMZYb+9QN8cHOyL1RuoKxoVLIRvzz6Cx+2/pCr8Vd5asNTvLfnPSJuW0HUm5RwfDnMbwZn1kOnSTBsB5RtrLeyAkUZAEWhIfxWOKO3jcbb2Zt5nefh6pC/p7iDV6JoVrEYdrnc+jGvONprk8t7L1pBR2hh7IQdA6oNYMOADQytM5T1Qet59PdH+fbkt6QYUh5cgSWIDdbSNf8xDIpXheGB0OFtbWVvEUMZADOT7lvLDqPRyM2b1psQzFZJMaTwxo43iE+J58vOX+Lnlr/c/REJyQSFJ9LUzNE/mWldpTgXwxO5GVc0NmHxcPLgzaZvsrrfapqXbs6nhz9lwJoB7AreVXAijEY4uFjbqOXybugxC/63GUrWLDgNVoYyAGbg4sWLrF69mhEjRlC2bFnWrl1751yTJk1IS0tj8ODBnDp1iuDgYJ588klu3rzJlClTmDZtGtOmTePIkSM6fgLbRkrJjP0zOB5+nGltpmWZ1jm3HLocDWD28M/MtKiszQMcNLVXVKjgVYF5XeaxsOtC7IQdo/4exYi/RnAp9pJlG464oG3NuOEtLUf/yH+g5QibS95mbpQBMAPTp09nz5497N+/nyNHjtCvXz+io6NJSUm5k/lzwoQJ7Ny5E1dXVz744APc3Nx49NFH6d27N3Z2dhw7dkzvj2GzrDy38s4+vo9UfOSh6jp0OQpnBzvqlrVs6F8dfy9cHO04fKVoGYB02pZty+99f2ds07Ecu3mMgWsGMvfgXOJT4s3bkCENdn8GX7WBGyeh33x4brWWt1+hFoKZg2+//ZakpCQOHjxI6dKlkVIybtw4XnrpJQDWrl2Lo6MjI0eOpH379rz66qv07dsXOzs7PDw8ePzxx1VK6Hxy+MZhZu6fSbuy7RjVcNRD13fwSjQNyvng7GDZJ0NHezvqB/hw+GrRNABwN63Eo5UfZd7Refx4+kfWBa3j9cav069qv3xN4N9D2AlYMwquH4eavbXkbZ4q425GCpUBmH1gNmeizpi1zpq+NRnXfFy258+cOcP48eO5desWp06d4rHHHqNYsWK4uLjcWdzVvHlznnnmGezt7XFycqJ///60bduW8ePHs2TJEn755Rezai4qhCWG8eaONwnwDGBW+1l3NnHPL7dS0jgVEssrHaqYSWHONKlQjG92BZGUajD7gjNbooRrCaa0nsITNZ5g5v6ZvL/3fZafXc745uNpWLJh3itMS4ZdH8HuT8G1GAz6AWr3KzT5e8yJcgE9JDVr1mT16tV07twZo9HIhx9+eN9uX6VLl2br1q088ojmnjhx4gQDBw5k06ZNvPPOO3rItnmS0pJ4ffvrJBuS+bzT52bJRnnsagxpRknTisXMoPDBNClfjDSj5Pi1mAJpz9qpU7wOS3suZWa7mYTfCue5P59jQuAEbt7KQ9DEtQPwVTvNANQbpCVvq9Nfdf7ZUKhGADk9qVuSK1eusGHDBipWrMioUaNYtGjRfWWOHTvGxYsXAahfvz6rVq1i/PjxzJgxg5UrVxa0ZJtn1oFZnIo8xRedvqCyj3ly0h+8HI0Q0LhCwRiA9HYOX42+Mylc1BFC0LtybzqX68ziE4v5/tT3/H31b4bVH8ZztZ/D2T6b1NzJCbBtGuz/CrwD4JlVUK1rwYq3QdQI4CGJjY3l2WefZe7cuXh6evLNN99QokQJkpKSMBjupsddvHgx5cuXB+DTTz/l8ccfZ/PmzXTr1o3YWLXfal5YfWH1nUnfTuU7ma3eQ1eiqFnaCy+XgokH93V3orKfO0eK6ERwTrg5uvFq41dZ038Nrcq04vMjn9N/dX+2X91+f5j1xW2wsBXsXwjNX9YifFTnnyuUAXhIgoKCePPNN2nQoAEA1apVw9fXl8WLF9OiRQuklISHh3PgwAFatGgBwBtvvMGePXsICwtjx44deHh46PkRbIqzUWeZvm86LUq3MMukbzpGo+TYtRgal/cxW525oUn5Yhy+Em31+wTrRTnPcnze+XO+7vY1TvZOvLr9VV756xWCYoLgdjSsHgVLB4C9MwzdBL0+AmdPvWXbDBYzAEKIJUKIvUKISdmcLyaE2CiECBRCfGUpHZamUaNGDBgwAIPBcF8kz4svvki5cuU4evQoQ4cOBSAlJeWeclu2bOHTTz+lTp06BarbFklISeCtnW/h6eRplknfjFyKTCQ+KY0G5XzMVmduaFKhGNG3UrkUkVig7doarf1b81vf3xjXbBwnwk/w2NoBzP6+NXEnfoW2b8Iru6FCK71l2hwWmQMQQgwE7KWUrYUQC4QQ1aSU5zMVew5YJqX8WQjxkxCiqZTykCX0FATu7u7s3r37nmNLliy58zp9AnjXrntXPvbo0YMePXpYXqCNI6Xk/b3vExwfzJLuSyjhat5EbekTsQ11MAAAh65EU9lPjQRzwtHOkWfLdaPXiQ3Mi/yHnzw92VilBmPK12GAvSNFN44q/1hqBNARWGF6vQ1om0WZSKCGEMIHKAdczaoiIcQwIcQhIcSh8PDwLBtTw+fcYcv36af/fmLrla283vh1mpRqYvb6j1+Lwd3JnioF3AlX8fPA08WBYyoSKGekhGO/wPzm+J77mw8ajmH5oz9TsVg1pvwzhcEbBnP0ptpVL69YygC4AyGm13FAVlsf7QaqAa8CZ4AsZ8KklIuklE2llE39/O7P7+Li4kJkZKRNd24FgZSSyMhIXFxc9JaSZ47dPMbHhz6mc7nODKkzxDJtXIuhXoA39hZKAJcddnaC+gHe/BscU6Dt2hQxV2HZY7D6FfCrqbl72r1FLb96fN/je+a0n0NUUhTP//k87+x6h7DEML0V2wyWCgNNANJTMXqQtaGZAbwipYwTQrwJDAXuj598AAEBAQQHB5Pd6EBxFxcXFwICAvSWkSeikqIYu3Mspd1LM7XtVIQF4rmT0wycvh7Hi23NE06aV+oH+KgFYVmRnrztr8na+54fQbOXIMPueUIIelbqSYeADnx78lu+O/kdO67t4KV6LzGkzpDsw0YVgOUMwGE0t88+oAFwNosybkA9IcQ+oAXwV34acnR0pFKlSvnVqbBiDEYDEwInEJ0UzbJey8yy2Csr/rseT6pB0rCct0XqfxANArxJM0r+ux5Ho/IFswbB6ok4D2tGw7V9UKUL9PkMfMpnW9zN0Y3RjUYzoNoAPj70MfOOzuP387/zdtO36Vy+s0UeHAoDlnIBrQaeE0J8AjwBnBJCTMtUZibaE38s4AuofAiKe1hycgl7Q/cyocUEahWvZbF20ieACzoCKJ36AVq7/war9SAYUiHwY1jYBsLPQP+F8OyqHDv/jJT1KMsnHT9h8SOLcXVw5fUdr/Py1pe5EH3BwsJtE4uMAExunY5AN2COlDIMOJ6pzAFAxT4qsuTIjSPMPzafXpV68Vi1xyza1vFrMZT0dKa0lz7zI2W8XSjh4czxoj4PcP249tQf9q+Wu6fnR+CZ1fThg2lRpgUr+6xkxdkVzD82n8fXPc6TNZ5kZMOReDvrM9KzRiy2DkBKGS2lXGHq/BWKXBOTFMM7u94hwCOA91u9b/Hh+7FrMTQo56Obm0AIQYMA76I7AkhNgr+mwKJOEB8GTyyFJ37Md+efjoOdA0/Xepr1A9bzePXH+fXsr/T+ozcrzq7AYDQ8uIIigFoJrLAqpJS8t/c9IpMimdNhDu6O7hZtL/ZWKkERiQUe/5+Z+gE+XAxPICG5iKUFv/IPfNUWdn8CDQbD6ANQu69ZmyjmUoxJLSexovcKqvpUZeq+qTy5/kkOhdnssiOzoQyAwqr4+czP7Li2g7eavEWd4pb3EP4bEgMU/AKwzNQv542UcKKojAKS42HDWPiuBxiS4bk/oP98LX2zhajhW4Nvu3/LRx0+IjYllqGbhzJ251iuJ1y3WJvWjjIACqvhdORpPj70MR0DOvJMrWcKpM30CeB6Afr6hRvcmQiO0VVHgXDhL1jQSgvxbPEKjPgHqnQukKaFEPSo2IO1/dcyssFIdlzbQd/VfVl4fCFJaUVjf+aMKAOgsAoSUxN5e+fb+Lr4MrWNZeL9s+LYtViq+LkXWAbQ7PB1dyKgmGvhnge4FQV/vKIt6nJ01TZk7zkbnAs+BYargysjGo5gbf+1tA9oz4JjC+i3uh9bLm8pUotKlQFQ6I6Ukqn7phKcEMzs9rPxcfEpsLZPhMTcCcPUmwYBPoU3EujUapjfHE6shHZjYXgglG+htyr8Pfz5uOPHfNv9WzycPHhr51u8uOVFzkWf01tagaAMgEJ3Vl9YzYagDYxsMNIieX6y42Z8Ejfikqlb1jrCAusHeBMcfZvIhGS9pZiP+DBY/iysHAJe/vDydujyHjhaV0qSZqWbsbz3cia1mMS56HMMWjeIafumEZMUo7c0i6IMgEJXgmKCmHlgJi1Kt+Clei8VaNunQuMAqOtvmRXGeeXOgrCQQuAGkhKOLtOe+s9vha5T4KVtUKa+3sqyxcHOgSdrPsmGARt4ssaT/HbuN3qv7s0vZ34hzVg4o7OUAVDoRlJaEmN3jcXVwZWZ7WaaNb9/bjhl6mhrW4kBqBfgjRDw7zUbNwDRl2Fpf1gzCkrWgVf2QNvXwd42dqD1dvbm3RbvsrLPSmoWq8mM/TN4Yv0THLh+QG9pZkcZAIVuzD00l/PR55nedjp+bvdnerU0J0PiqFTCHU+dJ4DT8XB2oHIJd06G2qgBMBpg31dahE/wIXj0Y3hhA5SoqreyfFGtWDW+eeQbPun4CYkpiby45UXe3PEmoQmhekszG7ZhkhWFjq1XtrL87HKG1h1K27JZbRdheU6Gxuoe/5+ZOv7eHLbFPYLDz2ppHIIPQNVu0PtT8Cmnt6qHRghBtwrdaFe2HT+c+oElJ5ewK3gXQ+sO5X91/4erg+uDK7Fi1AhAUeCEJITwwZ4PqF+iPmMajdFFQ8ytFIKjb1vNBHA6dct6ERJzm6jEFL2l5A5DKuz6SFvNG3keBiyCZ1YWis4/Iy4OLgxvMJy1/dfSuVxnvjr+FX1X92XTpU02HTaqDICiQEkzpjEhcAISyez2s3G008f9cncC2LoMQB2TnlO24AYKPQqLOsK2aVDzURh1EBo8CYU49XJp99LM6TCH73t8j4+zD2/vepuhm4dyJuqM3tLyhTIAigLlmxPfcPTmUSa1nESAp36b05w0TQDXsZIJ4HTS9aQbKKsk9TZsfR++6QyJEfDkTzDoe/Ao+HkcvWhSqgm/Pvor77d6n6CYIJ5c/yQf/vMh0Um25b5TBkBRYBy7eYyvj39N78q9ebTyo7pqORkaR1kfV4q5O+mqIzM+btqK4JPWGgp6eY+Wq3/P59DoWRi1H2r11luVLtjb2TOo+iDWDVjH0zWf5vfzv/PoH4/y038/2UzYqDIAigIhPiWe8YHjKe1emoktJuoth1MhsdQta11P/+nU8ffitLWNAJLiYP2b8H0vMKbB82ug7zxw9dFbme54O3szrvk4VvVdRZ3idZh1YBaD1g1i3/V9ekt7IMoAKAqE6funE5YYxuz2s/FwKvjcLxmJT9JSQFub/z+duv7eBEUkWk9q6HNbtNDOQ99Cy1Ew8h+o3FFvVVZHFZ8qLOq2iM86fcbttNu8vOVlXt/+OsHxwXpLyxZlABQWZ93FdWwI2sCIBiNo4NdAbzn8dz0egDrWOgIw6frvus6jgMRI+H0Y/DxIS9j24lboMQOcLLtHgy0jhKBL+S6s6b+GVxu9yt7QvfRb3Y8vjnzBrdRbesu7D2UAFBblWvw1pu+fTuOSjQs81UN2pPvXrXkEAOg3DyAlnFylpXE4uQo6jIPhu6BcM3302CDO9s68XP9l1vVfR7eK3fjmxDf0Wd2HDUEbrCpsVBkAhcVID/m0w45Z7WYVeKqH7DgZGoufpzMlddoD+EGU9NL2CNYlEijuOvz6NPz2Py2Wf/gu6PQuODgXvJZCQCn3UsxqN4ulPZdSwrUE4wPHM2TTEE5HntZbGqAMgMKCfP3v1xwPP877rd+njEcZveXc4VRInNUkgMuOumW9CnYEICUc/gHmt4CL2+CRafDiX1DK8ruyFQUalmzIL4/+wpTWU7gSd4Wn1j/F5L2TibwdqasuZQAUFuHwjcMs+ncR/ar0o0fFHnrLucPtFAPnb8Zb3QrgzNTx9+L8zQSSUgtg8/KoS/BjX1j3KpSuByP2QusxNpO8zVawE3YMrDaQ9QPW81zt51hzYQ19/ujD0tNLSTWm6qNJl1YVhZq4lDgmBE6grEdZJrSYoLecezgTFodR3l1xa63U9ffGYJScuxFvuUaMBvhnvhbhE3IUen8GQ9ZB8SqWa1OBp5Mnbzd7m1X9VlG/ZH3mHJzD42sfZ2/I3gLXogyAwqxIKZn6z1TCb4Uzu91s3B2tK2LkZHoKCCuNAEqnzp2JYAvNA9w4DUsegc3vQqX22oKupkPBTnUJBUVl78os7LKQLzt/SaoxleF/DWfMtjFci7tWYBrUX1thVtZeXMumy5sY1WgU9fzq6S3nPk6HxuLj5khZH+vO4ljO1xVPFwfz5wRKS4Eds+Dr9hB9CR5bAk8vB++y5m1HkSuEEHQo14HV/VbzeuPXOXD9AP3W9OOzw58VSNioMgAKs3E17ioz9s+gaammDK0zVG85WXIyJI66/t4Ftul8fhFCUMff686IxSyEHIZFHWDHTKjTH0YdgHqPF+rkbbaCk70TL9Z7kfUD1tOzUk+WnFxC7z96s+7iOozSaLF2lQFQmIVUYyrjdo3Dwc5Bl929ckNKmpGzYfFWuwAsM3X9vTlzPY40w0N2ACm3YPNEWNwVbsfA4F/hscXgXsIsOhXmw8/Nj+ltp7Os1zJKuZXi3d3v8tyfz3Ey4qRF2lMGQGEWFh5byMnIk0xuPZnS7qX1lpMl52/Gk2IwWu0CsMzUKetFcpqRi+GJ+a/kUiAsbA3/fAmNh8CofVCjp/lEKixCA78G/PToT0xtM5WQ+BCe3fisRXYiU3FeiofmYNhBFp9YzMBqA+lWoZvecrLllGlC1dpSQGdH3Qx7A9Qo7Zm3i5NitZTNh7+HYpW06J5K7c0vUmEx7IQd/av2p2v5ruwJ3YO/h7/Z21AGQPFQxCbHMiFwAhW8KjCu2Ti95eTIydBYPJwdqFjcuiKTsqOynwcujnacDIljYOM8XHj2Ty1zZ0KYFs/f8V1wcrOYToVl8XDyoHvF7hapWxkARb6RUjLlnylEJkWyrNcy3Bytu5M5FRpH7TJe2NnZxqSnvZ2gVhmv3G8SnxgBf46Dk79ByTrw1DIo28SyIhU2jZoDUOSb1RdWs/XKVl5t9Cp1ilt3ygCDUXI6NI7aNuL+Saeuvzf/hcZhNOaQQExK+HclfNkMTq/RnviH7VCdv+KBWMwACCGWCCH2CiEmPaDcAiFEH0vpUFiGy7GXmXlgJi3KtGBInSF6y3kglyISuZ1qsPoUEJmp4+9FfHIaV6OyiQmPDYFfnoLfXwLfyvBKIHQcBw7WtdOZwjqxiAEQQgwE7KWUrQF/IUS1bMq1A0pLKddZQofCMqQaUhkXOA4neyemt5mOnbD+gWT6gipbmQBOJ91g3ZcZ1GjUNmiZ3wKCdkL3GfDiFihZSweVClvFUv+5HYEVptfbgLaZCwghHIFvgMtCiH7ZVSSEGCaEOCSEOBQeHm4JrYo8Mu/YPE5HnmZK6ymUci+lt5xccSo0DicHO6qW1Hc3srxSrZQHjvbi3nmAyIta8rb1b0DZRtoOXa1GgRWuvVBYN5YyAO5AiOl1HJBVL/E8cBqYAzQXQozJqiIp5SIpZVMpZVM/Pz+LiFXknv3X9/P9ye8ZVH0QXcp30VtOrjkVGkvN0p442lv/aCUjzg72VCvpqaWGNqTBni+0uP7rx6HPF/D8WvCtpLdMhY1iqf+GBCA92YpHNu00AhZJKcOAZUAnC2lRmImYpBjeDXyXit4VebvZ23rLyTVSSk6GxNmc+yedumW9SAn5F7mkK2x9D6p01pK3NRmi0jgoHgpLhYEeRnP77AMaAGezKHMBqGx63RS4YiEtCjMgpeSDvR8QnRzNl12+xNXBupOpZSQk5jaxt1OtPgV0lqQl88ztn6htWIwx2gf7x7+DOgNUx68wC5YyAKuBQCGEP9ATeEoIMU1KmTEiaAnwrRDiKcAReNxCWhRm4Lfzv7Ht2jbGNh1LreK2NdF40sZWAN/h2kFYO5oG4Wf43dgW30c+pmPdmnqrUhQiLGIApJRxQoiOQDdgjsnNczxTmXhgkCXaV5iXoNgg5hyYQ2v/1jxX+zm95eSZ06GxdxZV2QQpibBtOuxbAF7+JD3xK28tNfJqhB0d9damKFRYbCWwlDKau5FAChslxZDCuF3jcHVwZVqbaTYR8pmZU6FxVPFzx8XRBqJkgnbA2lch5go0fRG6TsbFxYsqfjv12SReUahRqSAUOfLFkS84E3WGLzt/iZ+bbUZhnQyNpXUVK099fDsGtkyCo0vBtwq8sBEqtrlzuo6/FwcuRemnT1EoeeDjnBDivsBpIYSzEGKqEKJyVtcoCgd7Q/fyw+kfeLLGk3Qo10FvOfkiPD6ZG3HJ1u3/P7NBW9B17Gdo8zqM2HNP5w9aSojrsUlEJCTro1FRKMnNeP49IcS8TMdcgJvAAvNLUlgDUUlRTNw9kSreVRjbdKzecvLN3RXAVhgBlHATVr4Avz4N7n7w8t/QbQo43h9hlb6JjXIDKczJAw2AlHIcECGEWJThWKyUch5g3ekfFflCSskHez4gLjmO2e1n4+LgorekfJPeYVpVEjgp4fhymN9ce/rvPAmGbQf/RtleUifD3gAKhbnIcQ5ACFEJMADfA+8LIaYB6YbAD7DcZpUK3Vh+djk7gncwvvl4avjW0FvOQ3EqNJbyvm54uzrqLUUj5pqWwuHCVghoDv2+BL8H32NvV0fK+bre2dRGoTAHD5oEno1mACQggEeAfsAJ07GPLKpOUeBciL7A3ENzaVu2LU/XfFpvOQ/NqdA46lrDHsBGIxxaAn9NBmmEHrOh+ct5yt9T198793sDKBS5IEcDIKV8IuN7IYQX8BfwqpQywpLCFAVPsiGZcYHjcHd0Z2qbqQgbX20al5TKlchbPNG0nL5CIi7A2jFwdS9U7gR9PoNiFfNcTd2y3vx5Moy4pFS8XKxkRKOwaXKcAxBCTDat5gW0BV7ANO6mcFAUIj47/Bnnos8xtc1USrhaedhkLjgdqvMKYEMa7P5US9528xT0WwDP/ZGvzh/uzmOcVhPBCjOR7QhAaI9/Z4EVJiMQAtxGcwUhhLAD3KSUrQpCqMKy7A7ZzbL/lvFMrWdoH1A4Ng8/dccA6BABFHYC1ozSsnbW7A2PfgyepR+qyvRN4k+GxNKycnFzqFQUcbI1AFJKCfwC/CKEeBp4HXhFSnmqgLQpCoiI2xFM3D2RasWq8UaTN/SWYzZOhcRS0tMZP0/ngms0NQl2fQR7PgNXX3jiR6id7XYXecLP05lSXs4qFFRhNnK1ElhK+bMQ4iigVqEUMqSUvL/nfRJTE1n8yGKc7Quws7Qw2gRwAT79X90Pa0dDxDlo8DR0nw5uvmZtoo6/twoFVZiN3KwErgcgpfxPShmU4fgISwpTFAw/n/mZwJBA3mr6FtWKZblzp02SlGrgQnhCwfj/kxNg4zvwbXdIvQ3ProIBC83e+QPU9ffiws0EbqcYzF63ouiRm5XAXwIIIbaafm8wHX/eUqIUBcO56HN8cugTOgR04KkaT+ktx6ycCYvHYJSW9/9f+BsWtIIDi7SwzpH/QNWuFmuuTllvjBL+C1NuIMXD86CFYP2A20KI8oBjpt+2lxZScYektCTG7RqHl7MXH7b50OZDPjNj8U3gb0fD5olw7CcoXg2G/gkVLB8Pkf55ToXG0bh8MYu3pyjc5BQF1B1tVy9HtAVh1Uy/a5p+ly0IgQrL8MnhT7gQc4Gvu36Nr4v5XRV6czIkDm9XRwKKWWDnstNrYeNYSIyAtm9Ch3HgWDDpMsr6uOLj5sipEDUPoHh4choBBAIHgWVSysFCiI2m37+bfu8qII0KM7MreBe/nPmF52s/T+uyrfWWYxFOhcZSx9/LvCOb+Btax//fWihdD55ZCWUamK/+XCCEUCuCFWYjJzdOfWAv0EwI8S1Qz/S7iel3ddNvhQ0RcTuCSbsnUdO3Jq81fk1vORYhOc3Af9fjqB/gY54KpdRSNc9vDuc2Q5f34eXtBd75p1OnrBdnw+JJSVOpuBQPR07rAPYJIWqh5f75ENgPzEJz/wDMAWxgiyVFOkZpZOLuidxOu83sdrNxsnfSW5JFOBsWT6pBUj/ADBPA0Vdg/etwcRuUawl954Ff9Yev9yGo4+9NqkFy/ma8daa5VtgMD8oFJIUQDsBoIBVoBERKKf8qCHEK87Ls9DL2hu7lvZbvUdmn8Gbz+DdYc4/Ue5g1AEYjHPwG/poCQkCvudoWjXb6xz7UNU0EnwyJVQZA8VA8KBfQy8Bg4JCU8h+0jd27CSGOCCHeLgiBCvNwJuoMnx35jE7lOjGo+iC95ViUE8GxFHN7iAng8HPwXU/48x0o31IL7Wz+slV0/gAVi7vj6ezACTURrHhIHrQSeA+wREppFEI8K6VcBowTQuxD7QVgM9xOu824XePwcfZhSusphS7kMzP/hsRSP8An75/TkAp7Poeds8HRDfp/BQ2e0kYAVoSdnaBuWe87Ix2FIr886JFmPnd3/fqfEMJOCLEMzSUUaVFlCrMx9+BcLsVeYka7GRRzKdyx47dTDJy7EZ93/3/oMfimE2ybCjV6wuiD0HCw1XX+6dQv581/1+NITlMrghX554G5gKSUCaaXRrS9gL8D/kHbEEZh5Wy7uo0V51YwtM5QWpZpqbcci3P6ehwGo8y9/z/1tvbEv+cLcC8BTyyF2n0tK9IMNAjwIdUgOXM9ngblfPSWo7BRHmQAMnbyFYFlptejAGchRGLmTWMU1sPNWzf5YO8H1PKtxZhGY/SWUyCcCI4ByF0I6JV/tORtkReg0bPwyDRwtY0RUvoI53hwjDIAinzzIBdQxvFvkJRyILAVGCalfBRIyPoyhd4YpZF3d79LsiGZ2e1n42hfNHaQ+jck9k7a5GxJjocNY+G7HmBIgedWQ7/5NtP5g7YiuLi7E8evqXkARf554CSwEMIFSAM8MlyzQgjxn5TyfxZVp8g3P576kf3X9zO51WQqeVfSW06B8W9wLA0CvLOfAD7/lxbXHxsMLUZA50ng7JF1WStGCEH9AG/+NY14FIr8kOMIQEo5CfAElgOvCCHWA0JK2RlYUAD6FPngdORpPj/6OV3Ld2VgtYF6yykwEpLTuBieQL2yPvefvBUFf7wCPz2mRfi8uAV6zrLJzj+dBuV8uBCeQEJymt5SFDbKg7KBbkCbB6gPTAdqAClCiI6AnRDCSUrZy9IiFbnnVuotxu0ah6+LL5NbTy70IZ8ZORUSi5TcGwEkJZxeDRvf1jJ4tn9b+3Gw/Y1vGgT4IKXaIlKRfx7kAhoI+AArgdVAf+AWMAG4jpYpVGFFzDk4hytxV1jSfQnezkVrlWj6wqg7u4DFh8GGt+DMeijTUNuQvXQ9/QSamXRD929wjDIAinzxoFQQyUKIBGCRaRHYEiFES+AZKeUstLkBhZWw9cpWVp1fxUv1XqJZ6WZ6yylwjgfH4u/tgp+HExxZquXrNyRDtw+h5Siwz9UOqDZDcQ9nyvq4clwtCFPkk9ysA0jkbvgnUsp9wL4HXSeEWALUAjZKKaflUK4UsElK2ShXihVZEpYYxuS9k6lbvC4jG47UW44uHL8WQ6dSt2FpfwjaARXaQJ8voERVvaVZjAbl1ESwIv9YJLmJEGIgYC+lbA34CyFy2mx2LmCBXTuKDgajgYm7J5JqTNVCPu2KnmcuIu4WXWJXMTn4RQg+DI9+AkPWF+rOH7T1DteibhOVmKK3FIUNYqnsVh2BFabX29B2FrsPIURnIBEIy64iIcQwIcQhIcSh8PBwc+ssFHx36jsOhB3g3RbvUt6rvN5yCp6bZ3D4oScfOC4lsUxLGLUPmllH5k5Lk3EeQKHIK5b6D3EHQkyv44BSmQsIIZyA94HxOVUkpVwkpWwqpWzq5+dndqG2zsmIk8w/Op/uFbvTr0o/veUULGkpsHMOfN0Op9hLvJk2Cpchq8A7QG9lBUa9st4IgVoQpsgXlpoVS+CuW8eDrA3NeGC+lDKmKIUqmpP0kM8SbiV4r+V7RSrkk5AjsHYM3DgJdQbyesQThBk8cXEqXBO9D8LTxZHqJT05cjVabykKG8RSI4DD3HX7NAAuZ1GmKzBKCLEDaCiEWGwhLYWWmQdmEpwQzKx2s4pOyGfqbdjyHizuArci4amfSRu4hN3XBY3L204qB3PSuEIxjl6NxmhU+RkVecNSBmA18JwQ4hPgCeCUEOKeSCApZXspZUcpZUfgmJTyJQtpKZRsvryZ1RdW81K9l2hSqonecgqGy7thYWvY+wU0eg5G7oOaj3L2Rjy3Ugw0Ku+jt0JdaFzeh7gkbRW0QpEXLDJellLGmVYLdwPmSCnD0HYTy658R0voKKxcT7jOlH+mUL9EfV5p8IrecixPUhz89QEc+haKVYTn10LlDndOH70aA0CjckVzBNCkgva5D1+JplopT53VKGwJi4VJSCmjpZQrTJ2/wkwYjAbGB47HKI3Maj+r8Id8ntsMC1rC4e+h1WgYsfeezh/gyNVoSng4Uc63aEYTVyrhTjE3Rw5fUfMAirxRtGbMCgFLTi7hyM0jzGg7g3Ke5fSWYzkSI2HTeDixAvxqwhM/QkDTLIseuxpDw3LFitYkeAaEEDSpUIzDaiJYkUcKf6B0IeJ4+HEWHFtAz0o96V25t95yLIOUcOI3mN8MTv0BHcbD8F3Zdv7RiSkERSTSuIJPweq0MhpXKEZQeCLRakGYIg+oEYCNkJCSwPhd4yntXrrwhnzGhWrJ285uBP/G0O9LKFUnx0uOXYsBiq7/P50mpgioo9ei6VzzvmU3CkWWKANgI8w8MJPQxFC+7/E9nk6FbKJPSjjygxbeaUjVtmZsORLs7B946dGr0dgJLSdOUaZ+gA8OdoLDV5QBUOQeZQBsgI1BG1l7cS0jG4ykUclCljMvKgjWvgqXA6FiO+jzORSvkuvLj1yNoWZpL9yK2AKwzLg62VPb30tNBCvyhJoDsHJCEkKYum8qjUo24uX6L+stx3wYDbD3S1jQGq4fh96faeGdeej80wxGjl6NvhMGWdRpXL4Yx6/Fkmow6i1FYSMoA2DFpBnTmBA4AYCZ7WbiYFdInnJvnIYl3WDLRC2kc+Q+aDo0z8nbTl+PIzHFQPNKvhYSals0qVCM26kGzlyP11uKwkYoJD1K4eSbf7/h6M2jzG43m7IeZfWW8/CkpcDuT2DXXHDxgseWQN3HIJ8T2gcuRQEoA2Di7oKwKOoFFO05EUXuUCMAK+XYzWN89e9X9Knch16VC8G2y8GHYVEH2DET6vSHUQeh3uP57vwB9l+KomJxN0p5uZhPpw3j7+NKGW8XDl5W8wCK3KFGAFZIfEo84wPH4+/uz7st3tVbzsORcgu2T4d9C8CjNAxeDjV6PHS1RqPk4OUoHqmtIl4y0qKSL7svRCClLJyhwgqzokYAVsi0fdMISwxjVvtZeDh56C0n/1zaBQtbwT9fQuMh2kYtZuj8AS6EJxBzK5VmFZX7JyMtKxcnIiGFi+GJektR2ABqBGBlrLu4jo2XNjK64Wga+DXQW07+SIrVYvqP/ADFKmlbM1ZqZ9Ym9pv8/y0qFTdrvbZOy8ra/dgXFEnVkjb88KAoENQIwIq4Fn+N6fun07hkY16qZ6PZsc/+CfNbwNGl0HqMlrzNzJ0/aBPApb1cimwCuOyoUNyNUl7O7AuK1FuKwgZQIwArIdWYyvjA8dhhx6x2s7DPxSpYqyIxAv58B06ugpJ14KmfoKxl9imQUnLgUiQtKhVXfu5MCCFoWbk4ey9GqnkAxQNRIwAr4evjX/Nv+L+83/p9yniU0VtO7pES/l0JXzaD02uh00QYtsNinT/Alchb3IhLppkK/8ySFpWKEx6fTFCEmgdQ5IwaAVgBh28c5psT39C/an96VDTPJGmBEBsM69+E85uhbFMteVvJWhZvds/FCADaVFH+/6xoWVkzjPuDoqjip+YBFNmjRgA6E5cSx4TACQR4BDCh+QS95eQOoxEOLoH5LbUcPt1nwotbCqTzB9hzIQJ/bxcqlXAvkPZsjUol3CnpqeYBFA9GjQB0RErJh/98SPitcJb2Woqbo5vekh5M5EUteduV3VCpg5a8zbdSgTVvMEr2XoykW61Syr+dDenzAP8EqXkARc6oEYCOrLm4hs2XNzOq0Sjqlqirt5ycMaTBns+1TdnDTkDfefD8mgLt/AFOh8YRcyuVNlVLFGi7tkbbqiUIj0/m7A2VF0iRPWoEoBNX464yY/8MmpVuxtA6Q/WWkzNhJ2HtaAg9CjUehUc/Bi99JqrT/f+tqyr/f060q64ZyMBzEdQs7aWzGoW1okYAOpBqTGXcrnE42jkyo+0M6w35TEuGbdO1HD6xwTDoey28U6fOHzT/f41SnpT0VPl/cqKMtyvVSnqw63y43lIUVowaAejAgmMLOBl5kk86fkJp99J6y8maawdgzWiIOAv1n4IeM8FN37DLpFQDBy5F8UyLCrrqsBXaVfPjp/1XSEo14OJopQ8ZCl1RI4AC5mDYQZacWMJj1R6jW4Vuesu5n5RE+HM8LHlEe/3MbzDwa907f4AjV6JJTjPSRrl/ckW76iVITjNy8HKU3lIUVooaARQgscmxjA8cTwWvCrzT7B295dzPxe2w7lWIuQrNXoIuH2h5+62EnefCcbQXtKisDEBuaFHJFyd7OwLPR9Cump/echRWiBoBFBBSSqb8M4WopChmtZ9lXSGft2NgzShY2h/sHOGFjdpErxV1/gDbztykeSVfPJzVc0tucHNyoGnFYuw6p+YBFFmjDEAB8ceFP9h6ZSuvNnqVOsXr6C3nLv+t15K3HfsF2r4BI/ZAxTZ6q7qPa1G3OH8zgc41Vf7/vNCumh9nwuK5GZektxSFFaIMQAFwKfYSsw7MokWZFgypM0RvORoJN2HFEFj+DLj7wct/Q9fJ4Gid2TW3nbkJQOeaJXVWYlt0qqm5fv423T+FIiPKAFiYVIOW5dPZ3pnpbaZjJ3S+5VLC8V9hfnM4uxE6vwfDtoN/I311PYC/z9ykUgl3lf4hj9Qo5UlAMVf+On1DbykKK0Q5Uy3MvGPzOB15ms87fU4pd53dFzHXYP3rcOEvCGiuJW/zq6GvplxwKyWNfUGRPNdShX/mFSEEXWuV4pcDV7mVkoabk/qXV9xFjQAsyL7r+/ju5Hc8Uf0JOpfvrJ8QoxEOfAMLWsKVf6DnHPjfJpvo/AH2XIgkJc2o3D/55JHapUhOM7L7fITeUhRWhnocsBAxSTFMDJxIZe/KjG02Vj8hEedh7Ri4+g9U7qQlbytmW0/S287cwMPZQe3/m0+aVfLF08WBradv8EgdK114qNAFixkAIcQSoBawUUo5LYvz3sCvJg0JwJNSyhRL6SlIpJR8sPcDopOjmd91Pq4OOkysGtJg7xewYxY4ukC/BdDwabCxzJBpBiNbTt2gQw0/nBzUgDU/ONrb0alGSbaduYnBKLG3s63vgMJyWOQ/SggxELCXUrYG/IUQ1bIo9gzwiZSyGxAG2NBOKDmz8txKtl3bxuuNX6emb82CF3D9X1jcGf6eAtUfgVEHodEzNtf5g7b3b2RiCo/Ws6Fd0qyQbrVLEZmYwrFr0XpLUVgRlhoBdARWmF5vA9oC5zMWkFIuyPDWD8gyTk0IMQwYBlC+fHlz6zQ7QTFBfHTwI1r7t+bZ2s8WbOOpSbBrDuz+DNyKwxM/Qu1+BavBzGw4cR1XR3s61VD+/4ehQw0/nOzt2HgijCYVlCtNoWGpMbU7EGJ6HQdkG/4ihGgFFJNS7svqvJRykZSyqZSyqZ+fdS9nTzGk8M6ud3B1cGVam2kFG/J5dR981RYCP4b6T8Ko/Tbf+RuMks2nwuhcsySuTiqZ2cPg5eJI++p+bPj3Okaj1FuOwkqwVA+VAKQ7vj2ya0cI4QvMA/5nIR0FyudHPuds9FmmtpmKn1sBGavkBNj4DnzbA9KS4NlVMGChVSRve1gOXIoiIiGFXsr9Yxb6NvQnLC5JJYdT3MFSBuAwmtsHoAFwOXMBIYQTmptogpTyioV0FBh7Q/by4+kfearGU3Qo16FgGr3wNyxoBQcWQfNhMPIfqNq1YNouADaeuI6Lo92d1ayKh6NrrZK4Otqz9nio3lIUVoKlDMBq4DkhxCfAE8ApIUTmSKAXgSbARCHEDiHEkxbSYnGikqKYuGciVX2q8lbTtyzf4K0oWD0Slg0EB2ctpr/XHHD2tHzbBUSawcifJzX3j1q8ZB7cnBzoWrsUG09cJ9Vg1FuOwgqwyH+WlDJOCNER6AbMkVKGAcczlVkILLRE+wWJlJL397xPXHIcX3X9ChcHC+9UdXoNbBgLtyKh3VvQ/h0tzLOQEXg+goiEZPo1LKu3lEJF3wb+rDseyp4LEXRUE+tFHos9Wkkpo7kbCVRoWX52OTuDdzK++Xhq+FpwZW18GGwcC/+tg9L1NV9/mfqWa09nVhy6RnF3J7X618y0r14CLxcH1h4LVQZAoVJBPAwXoi8w99Bc2pZty9M1n7ZMI1LC0Z+05G3ntmibtLy8rVB3/lGJKfz13w36NyqLo736ipoTZwd7ejfwZ+PJ68TeTtVbjkJn1H9XPkk2JPNO4Du4O7oztc1UhCUWWUVfgaUDYM1IKFlby9Xf7k2wdzR/W1bEmmMhpBokg5oG6C2lUDK4WXmSUo2sORby4MKKQo0yAPnks8OfcT76PNPaTKOEawnzVm40wv6vtQif4IPQa662S1eJrBZUFz5WHgqmXllvapa2rh3JCgv1ArypW9aLn/dfRUq1JqAoowxAPggMDmTZf8t4ttaztAtoZ97Kw8/Cdz3gz3egQisttLP5y2BXNP5UJ0NiOX09Tj39W5jBzctzJiye48GxektR6EjR6FXMSMTtCCbtmUT1YtV5vcnr5qvYkAq75mqreSPOwYCv4ZnfwMf601+Yk+/3XsbNyZ5+DVT0jyXp28AfNyd7ftl/VW8pCh1RBiAPSCl5b897JKYmMrvdbJztnc1TcegxWNQJtk2FGr1g1AFo8JRNJm97GMLjk1l7LJTHmwTg7Va45zn0xtPFkT71/Vl7PJTYW2oyuKiiDEAe+PnMz+wO2c3YpmOpWqzqw1eYehu2fgDfdIbEm/DkMnjiB/AomuF5P+2/QorByAutK+otpUjwQpuK3E41sGy/zS/EV+QTZQByydmos3xy6BM6BnTkyRpmWLR8Za/m7tnzmZanf9R+qNXn4eu1UZJSDSzbd4XONUtS2c9DbzlFglplvOhQ3Y/v9lwiKdWgtxyFDigDkAuS0pIYt2scXs5eTGkz5eFCPpPjYcNb8F1PMKTAc6u1vXldi5lNry2y9ngoEQkp/K9NJb2lFCmGd6hMREIKq44E6y1FoQPKAOSCjw99zMXYi0xvMx1fl4fIsnl+K8xvCQeXQMuRMHIfVOlkPqE2SprByFc7LlKrjBdtqhbXW06RolXl4tQP8OabXUEYVJroIocyAA9g57Wd/Hr2V56v/Tyty7bOXyW3ouD34fDT4+DkDi9ugR4ztdcK1v0bSlBEIq91qWaZBXWKbBFCMLx9FS5H3mL9vypLaFFDGYAcCL8Vznt73qOmb01ea/xa3iuQEk7+Dl82g5O/aYnbXgmEcs3NL9ZGSTMYmff3BWqV8eKR2tnuG6SwID3qlqZmaU8+3nKOlDSVJbQooQxANhilkUl7JnE77Taz283Gyd4pbxXEXYflz8JvQ8E7AIbtgM4TtfTNijv8evAaQRGJvN61GnZqs3JdsLcTjOtZk6tRt/hJRQQVKZQByIZlp5exN3Qvbzd7m8o+lXN/oZRw5EeY3wIu/AXdPoSX/obS9Swn1kaJT0rl063naF7RVz3960zH6n60qlycedsuEJ+k1gUUFZQByIIzUWf47MhndC7XmUHVB+X+wqhL8GM/WDsGSteFEXuhzWtgrzY0yYovt18gMjGFSb1rKd+/zgghmNCrJlGJKSzYcVFvOYoCQhmATNxOu807u96hmHMxprTOZcin0QD/LICFrSHkCDz6CQxZD8WrWF6wjXImLI4lgZd4vEkA9QN89JajAOoH+DCwcVkWBwZx7ka83nIUBYAyAJn46OBHXI69zPR20/Fx8XnwBTf/gyWPwOYJULEtjNoHzV4sMsnb8oPRKHn39xN4ujjwbq9aestRZGBir1p4ODvw7u8nMKqw0EKP6qUy8PfVv1l5biUv1H2BlmVa5lw4LQV2zoGv2kFUEAz8Bp5eoU34KnJkUWAQR67G8F7v2vi653FyXWFRins4M/HR2hy6Es2S3Zf0lqOwMMo5beJG4g0+2PsBtYvXZkzDMTkXDjkMa8bAzVNQ9zHoMRs8/ApGqI1zMiSWj7ecpWfd0gxopDJ+WiOPNS7L1tNhfLT5LK2rFqeOv7fekhQWQo0A0EI+J+6ZSIohhVntZuGY3Y5bKbdgy3uwuCvcjoKnfoHHv1Wdfy6JuZXCyJ+O4OvuxIwB9dTEr5UihGDWwPoUc3dk1E9H1NaRhRhlAIAfTv3A/uv7Gd98PJW8s8lFcykQvmoDe7+ARs9pydtq9ipYoTZMmsHImF+Ocj32NgueaUIx5fqxaoq5OzH/6caExNzm1V+OkmZQC8QKI0XeAJyKPMUXR7+gW4VuDKg64P4CSbGw7nX4oTdIIzy/Fvp+AS5qWJxbpJRM/OMkgecjmNqvLk0qFO3Ed7ZC04q+TOlbl53nwhmvJoULJUV6DuBW6i3G7xpPcZfifNDqg/tdEuc2a51/Qhi0Gg2dJoKTmy5abRUpJdM3/MfyQ9d4tXNVnmpetHY4s3WeblGem/FJfPbXedyd7Jnct45y3RUiirQBmHNwDlfirrCk+xK8nTM80SdGwKbxcGIl+NWCJ5dCQFP9hNooBqNk8tpTLN13hSGtKvBGt+p6S1Lkg9e6VCMxOY1vAi9hkJLJfergYF/knQeFgiJrALZe2cqq86t4qd5LNCvdTDsoJZxcpW3InhQHHcZDu7fAQfmr80rsrVReW36UHWfDGd6hMuN71FRPjjaKEIJ3e9XCzk7w9c4ggqNvM29wIzxd1Ladtk6RNABhiWFM3juZusXrMrLhSO1gXCisfxPO/Qn+jbVNWkrV0VeojXIyJJZRPx8hNOY20wfU5enm5VXnb+MIIZjQsxYVfN15b81J+s3fwydPNKRhOR+9pSkegiJnAAxGA+/ufpdUYyqz28/GEXs49B1sfR8MqfDIdGg5Auzs9ZZqc9xKSeOzv86zZPclSng48euwVmrCt5DxdIvyVCrhzlsrjvHYwr2M6FCFkZ2q4OZU5LqSQkGR+6t9d+o7DoYdZGqbqZRPTYUf+8LlQKjYTovu8c1D5k8FAMlpBlYcvMb87RcJi0ticPNyjO9RC2835SIojLSqUpxNb7Tnw3Wn+XL7BZYfusbrXasxqEk5nBzU3IAtIaS0ndCupk2bykOHDuX7+pMRJ3lu43N0Ld+FOY7lEdtngL0jPDIVGg8B5abIE9eibrHi0DVWHLrGjbhkmlUsxtvda9K80kNsm6mwKQ5djmLmn2c4fCWaEh7OPN2iPE82K0dZH1e9pSkyIIQ4LKW8L5KlyBiAxNREnlj3BKmpt/gt1ohX6DGo3hN6fwJe/uYVWkhJMxg5ERLL9rPh7Dh7k3+DY7ET0L66Hy+3q0zrKsWVr78IIqUk8HwE3++9zPazN5ESGpTz4ZHapWhZ2Zd6ZX3UyEBnsjMAFnMBCSGWALWAjVLKafktYy5m7ptOcPw1vg0Lx0u4wmNLtDw+qsO6j6RUA8HRt7gadYtrUbe5GJ7Av8Gx/Hc9juQ0I3YCGpUvxtvdazCgUVn81dNekUYIQfvqfrSv7se1qFus//c6G09c56PNZwFwcbSjXllvqpfypFpJDyr7eeDv40IpLxcVSaQzFjEAQoiBgL2UsrUQYoEQopqU8nxey5iLrzZ/xJqwdQyPjqVUiS5srjeBZGMx+Pc6WY2A0g9J5P3HMhSXd87J+45xTzmZw7XZl7tHmemgvP9Qlu1nrNcoJclpRpJTDdrvNCNJd14biLudRvStFGJupRJzK4XEFEPGlnF3sqdOWW+ea1mBBuV8aFu1hErloMiScr5ujOhYhREdqxCZkMzBy9EcuBTFv8ExrDseSlxS2j3lPZwdKO7hhKeLA57Ojni6OODh4oCnswNODnY42NvhaG+Hk73A0V5772QvsLezQwgQYPqtvdHei7vHTefSn/MynrM1utcpjaOZ119YxAUkhPgC2CSl3CiEeBzwlFJ+l9cypnLDgGEA5cuXb3LlSt72LE1LTaLXD03xMRhxu/IUO4xN8vuxCgUujnY4O9jj7GCHs6MdTvZ2eLo4UszNkWJuTvi4OVHcw4mAYq4EFHOjnK8rfh7OyrWjeGiklITHJxMUkciNuCTCYpO4HptEVGIK8UmpxCelkZCcdud3qsFo+rEdN7UlOTWlO+7O+XtmL2gXkDsQYnodB1TNZxmklIuARaDNAeRViIOjC3NbfEaic3FKeFdl0j39mPYmY9+W/jK9w8tY/M5TRIajWfWLGZ827q83Y1vZt08O5e7V9GCdwg6cHbTOXnXkCr0QQlDSy4WSXi55uk5KSZpRasYgTZJiMGIwSiQSKbVRr5TynlH6fedIHxVLbGja8x5cHc0fmm4pA5AApDuGPcg66VxuypiF+vW6WqpqhUJhYYQQOJpcQCjPo1mxVKd7GGhret0AuJzPMgqFQqGwEJYaAawGAoUQ/kBP4CkhxDQp5aQcyjxgD0aFQqFQmBOLjACklHFAR2Af0ElKeTxT559VmVhLaFEoFApF1lhsHYCUMhpY8bBlFAqFQmEZ1PI8hUKhKKIoA6BQKBRFFGUAFAqFoohiU8nghBDhQN6WAt+lBBBhRjnmQunKG0pX3rFWbUpX3ngYXRWklH6ZD9qUAXgYhBCHsloKrTdKV95QuvKOtWpTuvKGJXQpF5BCoVAUUZQBUCgUiiJKUTIAi/QWkA1KV95QuvKOtWpTuvKG2XUVmTkAhUKhUNxLURoBKBQKhSIDygAoFIoijRDCVwjRTQhRQm8tBY0yADohhHAQQlwVQuww/dTTW5O1IoQoJYQINL0uK4QIznDf7ottLsoIIbyFEH8KIbYKIf4QQjhZy/fMGjtaIUQZYAPQHNguhPCzlvtVEBSJOYCC3Hw+twghGgNPSinH6a0lHSFEKeA3KWU7IYQj8AfgCyyWUn6rk6ZiwC9ASSllY9Ne0qWklAv10JNBlzfwK1pCxQTgSWAhOn/PhBAjgfNSyq1CiIXAdcBd7++ZqaP9HVgPPAV0Bmah//3qCiRIKfcJIeYC4YCv3vcrHdP/5CYpZSNL9GOFfgSQcfN5wF8IUU1vTSZaAgOEELuFED8JISyWmTU3mDraH9C26gQYAxwy3bfeQghPnaQZ0DrXONP7lsBIIcQ/QohPddIE8AzwiZSyGxCG1qnp/j2TUi6QUm41vfUD0rCO71kd4A0p5XRgM5oBsIb79Zep82+PNgq4jXXcr3TmAq6W6scKvQFA23MgPeX0Nu7uQqY3B4EOUsq2QAzQS18593W0Hbl73/YCuqyMlFLGZdor4k+gtZSyFVBdCFFfJ12ZO9pnsaLvmRCiFVAM2IoVfM+y6Gi7YyX3S2gbZT8JpALHsYL7ZdLVGUhEe8DoiAXuV1EwAJk3ny+lo5aM/CulvG56fQbQdWSSRUdrrfdtr5Qy3vRa9/uWoaO9hpXcLyGELzAP+B9W9D3L1NEKrOR+SY1RaA86pa3hfgkhnID3gfGmQxb5fywKBqDANp/PI0uFEA2EEPbAALQnD2vCWu/bZiFEGSGEG9pT5Em9hGTqaK3ifpk6jhXABCnlFazoe5apo22JddyvcUKI501vfYCvrOR+jQfmSyljTO8t8v2yln9qS2Ktm89/CCwFjgH/SCn/0lfOfVjrfZsCbEfbSvQrKeVZPURk0dFay/16EWgCTBRC7ABOYQXfsyw62llYx/1aBDwnhNgF2APtsYL7BXQFRpn+hg2BPljgfhX6KCAhhBcQCPyNafN5tf9w9gghdkgpOwohKgAbgb+A1mj3zaCvOutBCDECmMHdJ8TvgDdR37MsMQUZrACc0UZtE4BdqPv1QExGoC8W6McKvQGAO1++bsAuKWWY3npsBSGEP9pTx2b1z/lg1Pcsb6j7lTcscb+KhAFQKBQKxf0UhTkAhUKhUGSBMgAKhUJRRFEGQFEoEUKUE0Lk6fsthPAXQvg8RJtCCNEyi+MVhBBV81uvQmEp9F7mrFBYilGApxDiE7QInVtoYX72gESLQKkFJAEuaFEpw4FNwI7MlQkhtkkpOwsh3gMeQ1sl6gOskVJ+YCrWFugP7BNCjEUL36uGlo9nDXDBVNeHaKs5uwHxwHzgN6BXxkgrIURZtPDSM5nk1ACaSSmD83FfFIo7KAOgKKxMAtYBRillewAhREego5RyshCiOdpq1GVoaRy80Dr2ltqCVQBeAoxARcBNCNENLYzxdSnlDiFEb7TOOJ0ngQ+EEC8B1U0aXgB2SCl3mDR4oK3kbA2URFvRWRFIlFIa0kctUkojkJLD51MhuYqHRhkARaFESpmGFi+dbRHT7x/Q8uUMAL6UUn4OIIRYgLby0h6ogNbxl0UzGpgWNT0PvG163waoazrWEi2rZFZ4A8WB0WiLjfaijVaqmhYjVUUbRRxA6+R/Bk5nqqMWORsHhSJXKAOgKHQIId5E69D9pJQ1symWbgCiTb8jgEghxCtSyq/Q/jfSpJQngGNCiMeAi0Bl4BPuJs37H1rm1FDgK7Sn/o7AZDTj4g88KoS4hWYcktGe+D8GagKlgfrARDQX0XAp5QEhxBDgaSDWVMc9HxH4SQjxuZTyzzzdHIUiA2oSWFHokFJ+IqVsB+TkIxeZ3u8EWgHHhRA90J7804QQjwghfkBz2XRHS6/8N9qT+WS0OQSAq2hpBJaZyglgCDAbmCGl7GBKGeGAls7CEc0IfIo2R9AEKA8EmT7DD8DraPMTh9AMQZzptStaamXV+SseCjUCUBRV7NDcOkbT+9porpXGQA8gCq2zl8C3QHEp5STT5O4BtM75HFouGdDmEEqjZQWN4u4I4z+gC9rmOqC5k6ahTQ43QJsoroS2WQpo6a4zamxiqrskmlGpbrrW/qE+vUKBGgEoCjlCiEZCiGeyOBWE1vmfQst5dAhtLuAi2sjBAUiVUm6VUu7MdG2SqXwzKWV6NtIYtEiicDTXjRvwKFqn3VUIUVoIUU9KuQctJ87faKOFn6SUqcARoB/aPhHpxKMZn9Wm44dMr7/jrgtKocg3ygAoCjNl0fz1O0xbXHpx98k8Ec3lswJYgBaFM0dKuUlKOQJtovZ2hrrSn7gd0EYOTYFUU1I4gBbAO2ibd4xAiw7qiZYs7k9gJVqGR9DmCv5Gmy+oK4SohLZjVgraCCSdzmjzCf2BZqY2+6O5mjKWUyjyhXIBKQolpiywR4BhUspEIcQ4YKjpB7Sn9LlSyqOmSePSQJgpDPNv4LCUMtJU13o0Hz9ovvt3gK/RNhP/SQgRJqX8Ay1FdXr7k4DdUspUIcRHaOsQvjUl9FqINgJpieZ6+g4YC9wAfhNCPI0WLfQcmu8ftBGFQMsLDzBGCBEipcw4YlAo8oRKBqdQPAAhhIMprDT9vZCmf5yMr/NTn2mXLLv0BWD5qU+hyC/KACgUCkURRc0BKBQKRRFFGQCFQqEooigDoFAoFEUUZQAUCoWiiPJ/OaG6iDPDEX4AAAAASUVORK5CYII=\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#4.2\n",
    "\n",
    "'''\n",
    "需求量x是随机的，衡量报童的收入，不应该是每天的收入，而是长期卖报的平均收入，即收入的数学期望\n",
    "设定每天需求量为x份的概率为p（x）（x=1,2,3,4...）则需要建立p(x)与a，b，c之间的优化模型\n",
    "设报童每天批发量n份报纸的平均收入为s（n），则有x<n,x=n,x>n\n",
    "分类讨论：\n",
    "    x<=n,则售出x份，退回n-x份\n",
    "    x>n,则售出n份，无退回\n",
    "最终的离散型概率约束函数为：s（n）=（0，n）sigma[(a-b)x-(b-c)(n-x)]p(x) + (n+1,oo)sigma(a-b)np(x) sigma为级数求和符号\n",
    "其中，p（x）以及a,b,c已知，求使s到达max的n\n",
    "可以对目标函数进行求导，确定目标函数与批发量之间的关系\n",
    "文章参考：https://www.docin.com/p-89200656.html\n",
    "视频参考：https://www.bilibili.com/video/BV1554y1e7Cj?spm_id_from=333.337.search-card.all.click\n",
    "'''\n",
    "\n",
    "'''\n",
    ":func 报童模型演算\n",
    "为便于计算，假设订购量满足x~N(20,3**2)的正态分布\n",
    "一份报纸售价为2元，一份报纸进货价为0.8元，如若没有售出回收价位0.3元\n",
    "'''\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']\n",
    "plt.rcParams['axes.unicode_minus'] = False\n",
    "\n",
    "def normal(x, u, v):#计算概率分布函数，也可以使用scipy中的模块\n",
    "    return 1 / (v * np.sqrt(2 * np.pi)) * np.exp(-(x - u) ** 2 / 2 / v**2)\n",
    "\n",
    "u = 20 #订购量均值\n",
    "v = 3#订购量标准差\n",
    "accu = 0.1\n",
    "x = np.arange(u - 1 * u, u + 1 * u, accu)#x属于0到40，即取0,40之间的数字间隔为0.1\n",
    "y = normal(x, u, v)#y为x所对应的概率分布函数\n",
    "\n",
    "price = 2#售价\n",
    "cost = 0.8#成本\n",
    "scrap_value = 0.3 #未售出的残值为0.3\n",
    "\n",
    "\"\"\"利润分情况考虑\n",
    "先考虑某一天的情况：\n",
    "    1. 进货量<=需求量则 profit = (2-0.8)*x*y\n",
    "    2. 进货量>需求量则 profit = 2*x-cost*x+0.3*（进货量-需求量）\n",
    "\"\"\"\n",
    "def profit(papersin, marketrequire):#单次利润计算（不考虑概率的情况下）\n",
    "    if papersin <= marketrequire:\n",
    "        return papersin * (price - cost)#此时的利润即为(2-0.8)*x\n",
    "    elif papersin > marketrequire:\n",
    "        return marketrequire* price - cost * papersin + (papersin-marketrequire)*scrap_value #即此时的利润为2*x-cost*x+0.3*（进货量-需求量）\n",
    "\n",
    "costline = x * cost # 订购x份报纸所产生的的成本\n",
    "\n",
    "profitline = [round(sum(np.array([profit(papersin, marketrequire) for marketrequire in x])*y),2) for papersin in x]\n",
    "'''\n",
    "该profitline列表推导式为整个编程的关键\n",
    "    外部的 for papersin in x 表示报童订购报纸的数量在不断变化从0到40不断变化\n",
    "    内部的 for marketrequire in x ，表示在订购报纸数量固定的前提下（相当于嵌套循环），市场需求量在0到40之间不断变化的情况\n",
    "        另外需要求平均收入，则需要乘以一个周期内概率分部函数y\n",
    "'''\n",
    "print(profitline)\n",
    "print(max(profitline))\n",
    "print(max(costline))\n",
    "\n",
    "plt.plot(x, y / max(y), label='正态分布') #概率密度函数如果不正则话，则无法有效可视化\n",
    "plt.plot(x, costline / max(costline), label='成本')\n",
    "plt.plot(x, profitline / max(profitline), label='利润')\n",
    "max_indx = np.argmax(profitline) #返回取得最大值时坐标的索引\n",
    "print(max_indx)\n",
    "print(max_indx, y[max_indx] / max(y))#输出最大坐标索引以及标准化后的数据\n",
    "show_max = '({x},{y})'.format(x=round(max_indx * accu, 2), y=round(profitline[max_indx] * accu, 2))#x轴和y轴同时保留两位小数\n",
    "print(show_max)#展示最大值时的横坐标以及纵坐标\n",
    "\n",
    "plt.annotate(show_max, xy=(max_indx * accu, profitline[max_indx] / max(profitline)))#标注文字说明，这两行就是说明和画点\n",
    "plt.scatter(max_indx * accu, profitline[max_indx] / max(profitline))#散点图，横轴标为最大值坐标索引*0.1,纵坐标为利润最大值时的标准化数据\n",
    "\n",
    "plt.title('订购报纸数量与利润,成本分布图')\n",
    "plt.xlabel('订购报纸数量')\n",
    "plt.ylabel('标准化')\n",
    "\n",
    "# 由于订货量为整数\n",
    "c = max_indx * accu\n",
    "print(c)\n",
    "for i, c in enumerate((int(c), int(c)+1)):\n",
    "    profit_max = sum(np.array([profit(c, marketrequire) for marketrequire in x])*y)*accu\n",
    "    print(c,profit_max)\n",
    "#所以应该选择22份报纸，较高利润为22.23元\n",
    "plt.legend()\n",
    "plt.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1175030974154046 \n",
      " 0.10369611951319053 \n",
      " 0.09151150428043267 \n",
      " 0.687289278790972\n",
      "平均一台家用电器收费为：2674.2934822234856\n"
     ]
    }
   ],
   "source": [
    "#4.3\n",
    "\n",
    "#求一台家用电器收费Y的数学期望\n",
    "import numpy as np\n",
    "from scipy import integrate\n",
    "var = lambda x: 1 / 8 * np.exp(-1 / 8 * x)\n",
    "v1,err1 = integrate.quad(var,0,1)\n",
    "v2,err2 = integrate.quad(var,1,2)\n",
    "v3,err3 = integrate.quad(var,2,3)\n",
    "v4,err4 = integrate.quad(var,3,10000)\n",
    "print(v1,'\\n',v2,\"\\n\",v3,'\\n',v4)\n",
    "E = v1*1500+2000*v2+2500*v3+3000*v4\n",
    "print('平均一台家用电器收费为：{}'.format(E))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     x   y          prop        t_days\n",
      "0    0  22  3.678794e-01  3.678794e+01\n",
      "1    1  37  3.678794e-01  3.678794e+01\n",
      "2    2  20  1.839397e-01  1.839397e+01\n",
      "3    3  13  6.131324e-02  6.131324e+00\n",
      "4    4   6  1.532831e-02  1.532831e+00\n",
      "5    5   2  3.065662e-03  3.065662e-01\n",
      "6    6   0  5.109437e-04  5.109437e-02\n",
      "7    7   0  7.299195e-05  7.299195e-03\n",
      "8    8   0  9.123994e-06  9.123994e-04\n",
      "9    9   0  1.013777e-06  1.013777e-04\n",
      "10  10   0  1.013777e-07  1.013777e-05\n",
      "11  11   0  9.216156e-09  9.216156e-07\n",
      "12  12   0  7.680130e-10  7.680130e-08\n",
      "Power_divergenceResult(statistic=36.21310241772298, pvalue=0.00015596555915410734)\n",
      "pvalue < 0.05 所以不满足均值为1的泊松分布\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\goodboy\\AppData\\Local\\Temp\\ipykernel_23712\\145443555.py:6: DeprecationWarning: Please use `Power_divergenceResult` from the `scipy.stats` namespace, the `scipy.stats.stats` namespace is deprecated.\n",
      "  from scipy.stats.stats import Power_divergenceResult\n"
     ]
    }
   ],
   "source": [
    "#4.4\n",
    "\n",
    "'''\n",
    "参考csdn文章链接：https://blog.csdn.net/weixin_45807161/article/details/121459735\n",
    "使用scipy.stats库中chisquare函数，需要数据量尽可能多，如果数据量过少则会出现报错\n",
    "'''\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from scipy.stats.stats import Power_divergenceResult\n",
    "d = {'x': range(0,13), 'y':[22, 37, 20, 13, 6, 2,0,0,0,0,0,0,0]}\n",
    "df = pd.DataFrame(d)\n",
    "Poiss=stats.poisson(mu=1)\n",
    "df['prop']=Poiss.pmf(df['x'])\n",
    "df['t_days']=100*df['prop']\n",
    "print(df)\n",
    "df1=pd.DataFrame(df)\n",
    "n=None\n",
    "for i in range(len(df1)):\n",
    "    if df1.iloc[i,3] < 5:\n",
    "        n =i\n",
    "        df1.iloc[i+1,:] = df1.iloc[i+1,:] + df1.iloc[i,:]\n",
    "    else:\n",
    "        break\n",
    "if n is not None:\n",
    "    df1 = df1.iloc[n+1:,:]\n",
    "result = stats.chisquare(df1['y'], df1['t_days'], ddof=1)\n",
    "print(result)#\n",
    "print(\"pvalue < 0.05 所以不满足均值为1的泊松分布\")\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "170.25\n",
      "[[172.]\n",
      " [169.]\n",
      " [169.]\n",
      " [171.]\n",
      " [167.]\n",
      " [171.]\n",
      " [168.]\n",
      " [165.]\n",
      " [169.]\n",
      " [168.]\n",
      " [166.]\n",
      " [168.]\n",
      " [164.]\n",
      " [170.]\n",
      " [165.]\n",
      " [160.]\n",
      " [175.]\n",
      " [173.]\n",
      " [172.]\n",
      " [168.]\n",
      " [155.]\n",
      " [176.]\n",
      " [172.]\n",
      " [169.]\n",
      " [176.]\n",
      " [173.]\n",
      " [168.]\n",
      " [169.]\n",
      " [167.]\n",
      " [170.]\n",
      " [166.]\n",
      " [161.]\n",
      " [173.]\n",
      " [175.]\n",
      " [158.]\n",
      " [170.]\n",
      " [169.]\n",
      " [173.]\n",
      " [164.]\n",
      " [165.]\n",
      " [167.]\n",
      " [171.]\n",
      " [166.]\n",
      " [166.]\n",
      " [172.]\n",
      " [173.]\n",
      " [178.]\n",
      " [163.]\n",
      " [169.]\n",
      " [169.]\n",
      " [178.]\n",
      " [177.]\n",
      " [170.]\n",
      " [167.]\n",
      " [169.]\n",
      " [173.]\n",
      " [170.]\n",
      " [160.]\n",
      " [179.]\n",
      " [172.]\n",
      " [163.]\n",
      " [173.]\n",
      " [165.]\n",
      " [176.]\n",
      " [162.]\n",
      " [165.]\n",
      " [172.]\n",
      " [177.]\n",
      " [182.]\n",
      " [175.]\n",
      " [170.]\n",
      " [170.]\n",
      " [169.]\n",
      " [186.]\n",
      " [174.]\n",
      " [163.]\n",
      " [172.]\n",
      " [176.]\n",
      " [166.]\n",
      " [167.]\n",
      " [172.]\n",
      " [177.]\n",
      " [177.]\n",
      " [169.]\n",
      " [166.]\n",
      " [182.]\n",
      " [176.]\n",
      " [172.]\n",
      " [173.]\n",
      " [174.]\n",
      " [171.]\n",
      " [175.]\n",
      " [165.]\n",
      " [169.]\n",
      " [168.]\n",
      " [177.]\n",
      " [184.]\n",
      " [166.]\n",
      " [171.]\n",
      " [170.]]\n"
     ]
    }
   ],
   "source": [
    "#4.5\n",
    "\n",
    "#利用卡方检验分析4.6身高数据是否符合正态分布\n",
    "import numpy as np\n",
    "a = np.loadtxt(r'C:\\Users\\goodboy\\Desktop\\12660014_Python数学实验与建模程序及数据\\04第4章  概率论与数理统计\\Pdata4_6_2.txt')\n",
    "h = a[:,::2]\n",
    "junzhi =h.mean()\n",
    "print(junzhi)\n",
    "h_1 = np.reshape(h,(-1,1))\n",
    "print(h_1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "卡方值:\n",
      " 12.093523375642192\n",
      "P值:\n",
      " 1.0\n",
      "自由度:\n",
      " 76\n",
      "判断变量:\n",
      " 0\n",
      "原数据数组同维度的理论预测值(预测结果):\n",
      " [[168.70343612 171.79160059 168.55400881 170.84522761 168.10572687]\n",
      " [167.310837   170.37350954 167.16264317 169.4349486  166.71806167]\n",
      " [165.71929515 168.75283407 165.57251101 167.82320117 165.13215859]\n",
      " [168.70343612 171.79160059 168.55400881 170.84522761 168.10572687]\n",
      " [168.70343612 171.79160059 168.55400881 170.84522761 168.10572687]\n",
      " [168.50449339 171.58901615 168.35524229 170.64375918 167.90748899]\n",
      " [165.71929515 168.75283407 165.57251101 167.82320117 165.13215859]\n",
      " [167.310837   170.37350954 167.16264317 169.4349486  166.71806167]\n",
      " [167.50977974 170.57609398 167.36140969 169.63641703 166.91629956]\n",
      " [169.49920705 172.60193833 169.34907489 171.65110132 168.89867841]\n",
      " [171.28969163 174.42519824 171.13797357 173.46431718 170.68281938]\n",
      " [169.89709251 173.0071072  169.74660793 172.05403818 169.29515419]\n",
      " [166.91295154 169.96834068 166.76511013 169.03201175 166.3215859 ]\n",
      " [173.27911894 176.45104258 173.12563877 175.47900147 172.66519824]\n",
      " [172.88123348 176.04587372 172.72810573 175.07606461 172.26872247]\n",
      " [167.9076652  170.98126285 167.75894273 170.03935389 167.31277533]\n",
      " [171.28969163 174.42519824 171.13797357 173.46431718 170.68281938]\n",
      " [174.47277533 177.66654919 174.31823789 176.68781204 173.85462555]\n",
      " [168.70343612 171.79160059 168.55400881 170.84522761 168.10572687]\n",
      " [172.68229075 175.84328928 172.52933921 174.87459618 172.07048458]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import chi2_contingency   # 列联表分析\n",
    "from scipy.stats import chi2               # 卡方分布\n",
    "def chi2_independence(alpha, data):\n",
    "    CMIN, p, dof, prediction_value = chi2_contingency(data)\n",
    "\n",
    "    if dof == 0:\n",
    "        print('自由度应该大于等于1')\n",
    "    elif dof == 1:\n",
    "        cv = chi2.isf(alpha * 0.5, dof)\n",
    "    else:\n",
    "        cv = chi2.isf(alpha * 0.5, dof-1)\n",
    "    if CMIN > cv:\n",
    "        re = 1    # 表示拒绝原假设\n",
    "    else:\n",
    "         re = 0   # 表示接受原假设\n",
    "    return CMIN, p, dof, re, prediction_value\n",
    "# 测试\n",
    "\n",
    "\n",
    "alpha1 = 0.05    # 置信度，常用0.01，0.05，用于确定拒绝域的临界值\n",
    "\n",
    "data1 = np.array([[43, 49,22,114], [8, 2,5,15],[47,44,30,121]])\n",
    "\n",
    "data2 = np.array([[43, 49,22], [8, 2,5],[47,44,30]]) # # 插入数据（用于测试）\n",
    "a = np.loadtxt(r'C:\\Users\\goodboy\\Desktop\\12660014_Python数学实验与建模程序及数据\\04第4章  概率论与数理统计\\Pdata4_6_2.txt')\n",
    "h = a[:,::2]\n",
    "h_1 = np.reshape(h,(-1,1))\n",
    "\n",
    "CMIN, p, freedom, re, prediction_value = chi2_independence(alpha1, h)\n",
    "\n",
    "print(\"卡方值:\\n\",CMIN)\n",
    "print(\"P值:\\n\",p)\n",
    "print(\"自由度:\\n\",freedom)\n",
    "print(\"判断变量:\\n\",re)\n",
    "print(\"原数据数组同维度的理论预测值(预测结果):\\n\",prediction_value)#肯定原假设，身高服从正态分布"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            df     sum_sq   mean_sq         F    PR(>F)\n",
      "C(x)       6.0   2.281309  0.380218  1.689811  0.138196\n",
      "Residual  63.0  14.175400  0.225006       NaN       NaN\n"
     ]
    }
   ],
   "source": [
    "#参考csdn文章链接：https://blog.csdn.net/weixin_45807161/article/details/121459735\n",
    "#4.6\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "y = np.array([1.13,1.07,4.04,4.07,4.05,4.04,4.02,4.06,4.10,4.04,\n",
    "              3.86,3.85,4.08,4.11,4.08,4.01,4.02,4.04,3.97,3.95,\n",
    "              4.00,4.02,4.01,4.01,4.04,3.99,4.03,3.97,3.98,3.98,\n",
    "              3.88,3.88,3.91,3.95,3.92,3.97,3.92,3.90,3.97,3.90,\n",
    "              4.02,3.95,4.02,3.89,3.91,4.01,3.89,3.89,3.99,4.00,\n",
    "              4.02,3.86,3.96,3.97,4.00,3.82,3.98,3.99,4.02,3.93,\n",
    "              4.00,4.02,4.03,4.04,4.10,3.81,3.91,3.96,4.05,4.06\n",
    "              ])\n",
    "x=np.hstack([np.full(10,1), np.full(10,2), np.full(10,3),\n",
    "             np.full(10,4), np.full(10,5), np.full(10,6), np.full(10,7)])\n",
    "d = {'x':x, 'y':y} #构造字典\n",
    "model = sm.formula.ols(\"y~C(x)\", d).fit() #构建模型\n",
    "anovat = sm.stats.anova_lm(model)\n",
    "print(anovat)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "拟合的多项式为:1.0091815032766387*x + -0.3877735575011099\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxLUlEQVR4nO3deVhV1frA8e+SUVFxwCk1adDKLLUw0zQcstRKvVl2rZtDpd7Spl9pZhoQOTReUssGb2qTpZaVZmneHFKbsHBMTcUJJxAQUZTp/f1x5MhmEgTOPhzez/OcB867z97nZXt8Way99lpGRFBKKVXxVbE7AaWUUmVDC7pSSnkILehKKeUhtKArpZSH0IKulFIewtuuNw4KCpLg4GC73l4ppSqk9evXJ4hIvYK22VbQg4ODiY6OtuvtlVKqQjLG7C1sm3a5KKWUh9CCrpRSHkILulJKeQgt6Eop5SG0oCullIfQgq6UUi4WHh5eLsc1ds22GBISIucbtnjmzBkSExM5ceIEWVlZLspMVQS+vr4EBQURGBhodypKlZgxhgutvcaY9SISUtA228ahn8+ZM2fYt28ftWvXJjg4GB8fH4wxdqel3ICIkJaWxoEDB/Dz88Pf39/ulJRyC27b5ZKYmEjt2rUJCgrC19dXi7lyMsZQrVo1goKCiI+PtzsdpYolPDwcY4yzluV8X5bdL27b5bJjxw6Cg4Px9fV1YVaqIsnIyGDPnj00b97c7lSUKpHy6nJx2xZ6VlYWPj4+dqeh3Ji3tzeZmZl2p6GU23Dbgg5oN4sqkn4+VEUVFhZWLsd164KulFKeqLyGLWpBr0BE5IL73ZRSnk8LuhvLysoiIyPD+fzo0aPUq1ePTZs25XvtyJEjmTBhQqHH2rVrF5MnTwbgoYce4qGHHkJEGDlyJIcPH2b79u1cffXVxMbGFpnToEGDWL9+/QX+RFbbtm3DGMMvv/xSJsdTqrLTgu5iU6ZMcQ5XOt/D29ubt956y7nvoUOHOHXqFFdeeaXlmKmpqcyZM4fatWsX+r5169Zl3rx5REVF4ePjg6+vL1999RWrV6+mfv36xMTEcPz4cZo1a1Zk/j/++CPHjx8H4NVXX6VWrVqFPn744Ycij+Xn52f5miM7O5tTp06RnZ1d5P5KKSst6C7m5+dH69atnd0nOY+RI0eybNkySywjI4ORI0c6942NjaVFixb5Rv9Mnz6d5s2b88QTTwCwaNGifO/7559/MnjwYAIDA4mLiyMuLo6NGzfSu3dvNmzYwMqVK+nTpw9Vqjg+EiLC6dOn8x3H29vb8pp+/fqRnJyc7xEcHIy3d9H3reVc1Lzuuussv8i8vLwICAhg48aNJTizSim3vVO0PIWHh5fbRYnzKWhkxunTp5k5cyaPPPKIJZ5TEHfu3MmIESM4fPgw8fHx3HHHHQCMHTuWli1bEhUVxeLFi/Hy8uKvv/5i4MCBvP7664wYMcJ5rEWLFrF161ZiY2PZsWMHl156KcYYMjMz+fnnn/nhhx/YtWsXM2bMcO7TunVrYmJiiI+PZ+7cufj5+ZGamsrixYvZvHkzSUlJJf5ZC7J69WquueYa5/OsrCxOnz5N/fr1i7W/UsqhUrbQIyIi7E4BgIEDB+Ln50fDhg3x9/fnpptucnZXeHt789xzzwFQo0YNevbsibe3NzfffDM9e/Zk9erVZGdnM3z4cB588EFCQhz3GVxxxRVMnTqVxx9/nJ9//tn5Xm+88QbTpk0jJSWFq666is6dO9O0aVMWLlxIu3bt2Lt3L6mpqYgI06ZNo0ePHvz222+A4xfOpk2b2LRpE+np6cTGxjq///TTTwvsbtmzZ0+++Xeys7MtrfiUlJQCz4uPjw+NGzfW+xCUKqFK2UJ3F35+fjzyyCNERUXl29alSxdn33KDBg0YPXo0s2bN4pFHHiE0NJSnnnqKzZs388UXX1C/fn0++OADTp06hYhQq1YtGjVqRP/+/YmJiaF+/fp8/fXXDBs2jHHjxrFz504yMzPZsWMHQ4cOpXbt2mRmZpKWlkZAQACHDh2ifv36zrt0mzZtyvvvvw/A4sWLeeyxx+jSpQtRUVHcd999zJ49G4D58+ezZ88eRo8eXeDPGxcXx8UXX5wvfvPNN1ue161bl4SEhAs9rUpVWpWmhe6KeRQuxIwZMwps4a5Zs8byumPHjvH333/Tvn17tmzZQs2aNXn44YeZMmUKH330EWvWrOHIkSOcOHGC/fv3s3PnTiZPnky9evVISUnh3Xff5c033+Shhx6iRo0atG7dmrlz59KhQwc++OADateuTUxMDABHjhzhkksuKTTn+Ph4hg0bxrFjxyzxffv2WUbgpKenWy5s5vyCWLFiRb5rCDmPWbNm6XQPSl2gStNCz91vXpp5FEorIyPDMqqjqBZ67mK4atUqrrnmGqpXr87atWtp164dvr6+PPzww5w4cQJwFOLcBgwYgDGGmjVrMn/+fHx8fJg3bx4ffPABO3fupEaNGlx//fX07duXmjVr8tNPP3HLLbfw119/0b1793w5bd++nZMnTzJixAjGjBlDcnIyc+bMYc6cOZbXffTRR87v//zzT9q0aQNQ7EKtXS1KXZhK00J3FwkJCc7hhadPn+bNN98scMjiqlWrSEtLc+63fPly6tWrR2ZmJt9++62z4I4bN45LLrmkwMfKlSud+zdu3Bg/Pz8eeOABjh49Ss2aNTHG0Lt3b+bMmcOtt97K119/zenTp/nzzz/p2LGjc98DBw5w/fXX88ADD5CRkcFnn33G2LFjOXXqFM8//7yzdf2f//yHwYMHO5+np6dz7bXXOo9T3GGIOj+LUhemRAXdGDPUGLMyT6ytMWadMeakMWa/MSbSGOPWvyjKax6F4ti+fTtNmzYF4LPPPnMWv5wbemJjY52xV1991bnfvffey8GDB7nhhhv44YcfuOeeewBHqzd3Ec19N2nuFvHBgwfZtm0bgYGBHD58GBHh3//+N/3798fX15fevXuzf/9+xo4dS5MmTSzj0Zs0acLMmTP57bffqFWrlvO4Bw8epGHDhjRv3pzg4GBefPFFFixYQHBwMMHBwSxZssQ5xBFwDoPs2rVroWPvhw4dypkzZ8rj1CvlFtbuW8uR1CPnf+EFKHbhNca0B97KE6sGfAn8AfQCJgH/BzxRhjmWObv6zUWE6OhoWrVqBTj6om+88UZ27txped306dN56qmnLN1CoaGhrF+/3llMf/31VwBLwczLy8vL+X21atWYMmUKNWrUYNasWXz55Zd8/PHHjB07FoCqVasyfPhw3nzzTf7973/nO1bbtm3zxX7//XeuvfZajh49ysqVK0lMTCQ1NZU9e/YQHBycr0V+/PhxEhISSEpKyvdITEx0fr97925tpSuPE3M4BhNh6DSrE5N+mlQu71Gsgm6M6QYsA7bl2dQBqAM8ISKrRWQG8AHQp0yz9BC//vorcXFxdOvWjTNnzjBw4ECys7Np1KiR5XWdO3fms88+Y8iQIZahfzt37iQmJoZJkybx8MMP88Ybb5RoxsH33nuPTz75hAMHDtC/f3/S09N5/fXXOXLkCFlZWRw8eBBwTBNwvmsMR44cISEhgZCQkEJzyB0/efIkN954I++++y4BAQEMHTqU2NhYatWqxbp167jrrrsIDAykVq1ajB8/3nmTlFIVXVpGGsFRwbR991yjaPzN48vlvYrbQu8MDAK+yRMPOnuM3P+j/QH9m7kAEydOpE2bNlx00UV0796dhIQEvvvuOxITE/njjz8ARxFs3bo1q1evZvny5URGRgJw4sQJBg4cyNChQ3nmmWf4+uuv6d69OyLCnDlz8nVdFCYuLo5vv/2Wp59+mo0bN1K9enX27t1Lly5d2LFjB0uXLuXDDz9k0KBBpKam5ts/MzOTgIAA5s+fzx133EG1atUKLf65W+izZ89GRBgxYgQ+Pj4kJiYydepUwHED0y+//MLcuXMB6NGjB++8845lHL1SFVHEygiqTarG3uN7Afj2vm+RMKFeQL3yecPCho/l6Y+tcvZrOLAyV7wJcAp4EagBdAVOAIPOd8zrr79eirJ169Yit1c0Bw8elJYtW8rnn38uWVlZ8tprr0lSUpKIiLz00ksCSPPmzSUjI8O5z9atW+XYsWMSHx8vHTp0kKZNm8qxY8csxx01apQMHDhQkpKSLA9Ali9f7nxdZGSkBAUFSWhoqPz000/OeFRUlPj7+8uoUaPk1KlTIiKybt06ady4sfTt29f5uk2bNsnLL78sw4YNk+PHj0vz5s1l3bp1IiLi6+srNWrUkMDAQOfDy8tLPv74YxEROX78uDRq1EhGjx7tPN6sWbPE399fUlJSRERk5MiR0qpVK+f22267Tdq2bSvZ2dlFnldP+5woz/Dbgd+EcJyPB7968Lyf5eICoqWwWl3YhgJfnKegn43dB0iux4wi9h8ORAPRF198cZFJe+J/1MzMzAL/URMTE2XHjh2F/oOfOnVK7rrrLvn111/zbRs2bJgMHjw4XxyQJUuWOJ/v2rVLNm3alO91x44dKzCekJCQ75dHYZKSkiQrK6vI1yxdulTi4+MLfd8DBw7I3r17nc83bdokMTEx531vT/ycqIor9Uyq1H+1vqWYHztVvP9HxVVUQS/RmqLGmHCgi4h0Ofs8CNgArAPmA1cCY4EXROS1oo51vjVF//rrL6666qpi51ZZnTp1CnBc9KyM9HOi3MVzy59jytopzuc/PPADt1x6S5m/T1Fripb2xqKhwFFgwNnfHBhj4oFXjDFvikhGkXurUqushVwpd7F231o6zerkfP5oyKO8dftbRexRfkpb0C8Htom1mb8JqI7jgumhUh5fKaXcUsqZFBq/0ZjUdMfgAX9vfw4/fZhA/0DbcirtDUDxQFtjjFeuWE8co1yKnltVKaUqqCe/f5LAKYHOYr56yGrSnk+ztZhD6VvoS4BxwApjzFqgOdAPmCki+VdHUEqpCmxF7Aq6fdjN+fzpDk/z2q1FXi50qVIVdBFZZ4zpi+MO0aeBk8CnQMHzpyqlVAWRe0K/pLQkgl4NIlsc91bUqVqHvU/upbpvdRszzK9EXS4iEp4zwiVXbJGIXCMiviJSW0QGiciJMs1SKaVcLGchnBGLRlDnlTrOYv7zQz9zbMwxtyvmUImmz/VUIlKi2//L+zhKeYzLwESc+z8xvvN4IrtF2pjQ+bn1rIieysvLiyVLllhiW7ZsYfr06c7nn376qWVe8YJkZmbSvn17Pv3002K/9/Lly+nYsWO+xSmmTp1Kr169SE9PL/axAJKSkko0kdb48eP5/vvvLbFZs2YxZMgQ0tPTycjQka7KPuHh4Zh6xlHIHzgbTIFxMs7tizloQbeFj48PAQEBgGNVn5SUFJKSknjyySfZvHkz4FjV57vvvivyOFFRUWzcuJGoqCjL3OlF2bx5MykpKdStW9cSX7FiBQEBAYUuQnH8+HEOHz5MYmKiZV3Qm2++mbvvvtsSS0pK4ujRo/nmgtm7dy9TpkzJF09MTGTZsmX4+vrStWtXXnvNfS4yqcpDRIgwETDqXGz98PXI68LE8In2JVYCWtBdJCsri+TkZLKysvDx8cHPz4+0tDTmzZvHFVdcQadOnejcuTPjxo0DHEW/atWqhR5v3bp1vPDCC3z//fc0b96ce++9t1gt5S1bttC3b19LLDk5maVLl3LdddexcuVK52PXrl3O10RFRdGsWTPnXOc5jy1btvDDDz/kizdr1oz//ve/lveZNGkS//znP7n77rstLfE6depQo0YNAF5++WUWLVqkc6Irlxr2zTCqvHiuHFbzqQbhcF2j6+xL6kIUNidAeT8q2+RcsbGxuee7EUBat24t8+bNkxYtWoiIyLJlywSQ6OhoefTRR+Whhx4q8FibNm2SoKAgeffdd0VEJC0tTUJDQ6VPnz6Smppa4D45E4DlfWzevFnefvtt8fPzk2bNmjkfAQEBlsm0CpKVlSVVq1aVpUuXnvfnX79+vdSrV08OHz4sGzZskIsuukjWr18vIo6Jui677DJJSUmRdevWyRtvvCG9e/eWP//887zH9bTPiXKtDYc3WOZdIRxJTksWEZGwsDB7kysERczloi10F2nWrBlpaWlkZWURGBhIdHQ0P//8s6Ul3qNHD9577z2uuOIKjDGWudBzrFmzhu7duxMWFsbw4cMB8Pf3Z/HixaSkpHDdddfx008/5dvP39+fm266ybmIxJ49ewDHqkbTp08nLCyMPXv2OB//+Mc/Cux+GTZsmHOKXi8vL9LS0rjtttssU/f26NHDss/+/fsZMGAAL730EgEBAURGRnL11VcTHx/P5MmTmTlzJrt376ZJkyaMGTOGnTt30qNHj3zdQkqVlWzJxkQYWr/T2hn75K5PkDBx3hxk9wLyF6LCjXJ58vsniTkcY2sObRq2IapnVIn2Mcbg7+/vfO7l5VVowQTHBc/co04yMzOZMmUKr732GpGRkQwYMIDk5GTLvh999BHPP/88oaGh9OzZk4cffpg+ffrg7e2Nl5cX3t7e1KpVy7LP+++/T2JiIo899pglnncx6xz+/v4MHjyY2bNnF/hzhoeHs379ektsw4YN7Nq1i1GjRjFmzBiysrLYunUrb7/9NomJiXTs2JGYmBiSk5MtP/PevXsLfA+lSmPgFwP5bPNnzucX1biIuP+LszGjslPhCrqnePfddzl06BBDhgwpcPuZM2ecBTU5OZnOnTvj7e3N6tWrefzxx3n88ccL3C8sLIw1a9bw3HPPMX36dPr161dkHgMGDKBv375ERETQq1cvunVz3AWXkZFh+QWUw8vLiwULFlgWoM4tOTmZLl26WGK9e/cmOjqaNm3acM899xAaGkrTpk2ZPHky4LhQ+9prr7F06VK2bt1KTEwMK1eu5OjRo5w8edKylJ5SF+q3uN9oP7O9JZb6XCoBvgE2ZVT2KlxBL2nL2F2cPHmSb775hjlz5pCSkkJycjIRERHOxaEBDh06RGJiIldffTWJiYnOtUdr1arF/Pnzueyyy/Dx8WHZsmV4e3vnW08050Kjj48Pq1atIi0tzfKaVatW5RtrHhQURHBwMLNmzeKll15yFvT09PQCW+gAd999d5Et9JiYGEusSpUqXH/99Xz55Zfs37+fUaMcwwhmzJjBvHnz2LRpEyLC2LFj6dChA+3bt2fEiBG0bt1ai7kqtazsLLwjraVu4b0L6XdlP3sSKkfah+4ix48fZ/To0XTr1o2AgACeffZZWrc+13934sQJ7rzzTmbOnAk4RqNcdtllzu1XXnklPj4+tGjRAj8/P7y8vPItO+fr64uvry/PPvssQL5RMgX1oecYNWoUK1asYMOGDYDjL4SCWujgGCNfq1atAh9TpkwpcJ8vv/ySIUOG8OSTT7JmzRoWL17M5ZdfzpNPPklMTAwdO3bkqaeeYsaMGcTGxjJlyhQd6aJKrfcnvS3FvGW9lkiYeGQxhwrYQq+oLrroIvbs2YO3tzeTJllX/D5z5gx9+vShdu3aTJo0iV27drF7926uvvrqfMfx9fV13ohTkC5duhTasi6oDz1HmzZtaNOmDXPnzqV169akpaVRvXrBtzbfd999JWqhA7z00kucPHmSsWPH0qRJE9q1a+dcUxSgQ4cOLFu2jGPHjjF16lSmTp1KnTp1CnwPpc5n9d7VhM4OtcTSnk/D37vgRoqn0ILuQt7e+U/37t272bZtGw0aNOCbb76hatWqTJgwAWMM48aNY9myZZZuh+J0QeTtiimuWbNmOX+JpKamUrNmzQs6TkEWLlxIgwYNLK3+06dPO5/ffffd3HTTTXzxxRfMmzfvvH3/ShUkIysD35esgw2+v/97brv8Npsyci3tcrFZUFAQN954I4sWLcLPz4+nnnqKxYsXs3LlSrZs2eKcIChHaeZbyczMdN7Nefz48Xzb27Rpg4+PD+C4ezMwsOC5nefMmZOvuyfnkTdfcPxyOHjwIAsXLiQ8PJw+ffrQoEEDZs2aBcCxY8d48cUXERFatWpFnz59APjf//6XsxatUud10wc3WYp5x6YdkTCpNMUctKC7VEJCAl999RUnT550Fs5BgwaxZs0aDh8+TLdu3fjwww/57rvv6Ny5M2+++SaTJ0+2dGGICEOHDi20oK5atYrs7Ox8752RkcHatWupXbs2tWvXJjg42BnP69ChQ8TFxdGkSRPA0ZL+66+/2L17NykpKfTv35/Y2NgCH0888QSnTp0iNjaW7du3c/LkSdLT0+nVqxejR4/m4MGD3Hvvvaxbt46hQ4fyzjvv0LJlS7Kysvjzzz+Ji4vj0Ucf5Y8//uDWW2/lm2++Kft/COVRlu1ahokwrNu/zhlLH5/O2gfX2piVTQq746i8H5XtTlERkXXr1km1atXkX//6l6SnpzvjX3zxhfj4+EinTp1k9+7dln1GjRolcXFxzuetWrWSWbNmFfoeoaGhMmbMmHzxyMhICQ0NdT4/fvy4NGvWTLZv3+6MRURESKdOnSQwMFAef/xxZzwmJka8vb0lICBAAgMDi/WoUaOG+Pj4yMqVK0VEZN++fZZ8zpw5I9dff73UqFFDoqKiJDs7W0RENmzYIA0bNhRArr32WsnIyCjijHrm50QVT1pGWr67PFfGrrQ7rXJHEXeKGrHpT9qQkBCJjo4udLunruaee3x5juzsbJYsWcLtt99+3i6VFi1aMG7cuEIvipbGli1b2LlzJx06dKB+/fplfvy8tm/fTp06dahXr54lnpCQwLx587jrrrto2LBhkcfw1M+JKlqrt1uxJX6L83nPy3vy3f1FT2bnKYwx60UkpMBtWtBVRaafk8pl/pb5DFgwwBLLnJCJV5XKc79CUQVdR7kopdxeypkUAqdYL9L/+vCv3ND4Bpsyck9uXdBFV9FRRbDrr0vlWrlXDQIIuSiE34f9blM27s1tR7n4+voWe9EGVTmlpaU5RwspzxP1S1S+Yn7zipu1mBfBbVvoQUFBHDhwgKCgIGrUqIG3t7e21hXgaJmnpaURFxdHgwYN7E5HlbHEtETqvmKdOnntg2vp2LQjJlxrQFHctqAHBgbi5+dHfHw8x44dK9G6lcrz+fj40KBBgzK9m1W5Rnh4eKFzjedrkTe7mVVDVrkgK8/gtqNclFKeyRiT7/pH5KpIXlj5giWW/UI2xhi6dOnCqlX5i3poaGih0zh7Mh3lopRyS4dTD9Po9UaW2B/D/6Bto7bO57mLdkG/DNQ5bntRVCnlOcLDw53TU4CjMJsIYynmd7a4EwkTSzFXJaNdLkoplzK3GrjJGsvpXjmfLl26VMpultyK6nLRFrpSyiX2Hd/nuOiZq5hveXQLElb8+00qezE/H+1DV0qVu7yjVx649gE+/MeHNmXjubSgK6XKzYhFI3jvj/csMQnTi5rlRQu6UqrM/X3sb1pMb2GJ7Xp8F5fWvtSmjCoHLehKqTKVt3tlVLtRTOs9zaZsKhct6EqpMvHPBf/k8y2fW2LaveJaJSroxpihwGAR6VLANn9gMzBfRJ4rm/SUUu5u89HNXDPjGkts/1P7aVKziU0ZVV7FLujGmPbAW8BvhbzkWSAAmFQGeSml3JyIUOVF68jncZ3GMbH7RJsyUsUq6MaYbsBCYFsh2y8BxgKjRORE2aWnlHJHPT/uydJdSy0x7V6xX3Fb6J2BQUBboEsB29/EUexnlU1aSil39Hvc79ww07pK0JFnjlA/oPzXoFXnV9yCHiki2caYfJMsGGN6AncCy4GPjTF/AW+IyMkyzFMpZaOCulcmd5/M2E5jbcpIFaRYt/6LSHYRmyef/doAaASEAWuNMVXzvtAYM9wYE22MiY6Pjy9xskop17tx5o35irmESZHFvLD5zlX5KtVcLsaYdkAb4H0RuVZEugLdgFbA4LyvF5H3RCRERELq1atXmrdWSpWz1XtXYyIMv8b96owlPZtUrL7yiIiI8kxNFaK0k3M1P/v1jZyAiKwG/gZal/LYSikXymlVZ2VnYSIMobNDndum95qOhAm1/GvZk5wqltIW9NSzX2PzxE8D6aU8tlLKhSIiImgxrQXekdZLaxImjLxh5Hn3L3DOc2O0+8WFSlvQowEhV2vcGFMfuALQpbmVslFJCunSnUshHP5O/NsZO/HciRINRQwPD0dEnCsK5XyvBd11SlXQReQgMBf4yBjzD2NMb+Br4AiwoAzyU0pdoOL0Y08In4CJMPT8pOe54EIIkzCq+1Yvx+xUeSiLuVweBCKBqUAdHK32XiJyugyOrZQqJ/VerUeCSbDEJEwc49RKKSysDA6iSqxELXQRCc87j4uInBGRMSLSVEQCRCRURAq8o1QpVb6K04/95V9fYiIMCafOFfO059MgnDKj3Sz20NkWlfIg4eHhzmJqjCH3msGnM09TdaL19pAF9yygf8v+gLaqPYEWdKUqgbxzlNfyr0XSs0mWmLaqKz4t6Ep5qLCwMD7a8BGDvhpkiaePT8fHy8emrFR50oKulAc6mX6SCBMBX52LLblvCb2a97ItJ1X+tKAr5WHydq9cVvsydj6+06ZslCtpQVfKQzy25DGm/z7dEsuckIlXFS+bMlKupgVdqQru6MmjNHitgSU2t/9c/tnqnzZlpOxS2lv/lVI2MhEmXzGXMCmwmOsoFs+nBV2pCuj+L+/P11ee9UJWkXOv6JS2nk+7XJSqQPYf38/FURdbYosGLuKOFnfYlJFyJ9pCV8pGJekGMREmXzGXMCmymOuUtpWLyX1rsCuFhIRIdHS0Le+tlLvIe3t+QXp/0pvvdn5niWW/kO0s0mX5Xsr9GWPWi0hIQdu0y0UpN/X3sb9pMb2FJbZi8Aq6BHexJyHl9rTLRSkXK043iIkwlmIe6BeIhEmpirlOvuX5tMtFKRvl7Qa54f0b+P2gdbGvkqwapDyfdrko5eY2HtlI63es66r/9vBvtGvczqaMVEWkBV0pG4WGhuYbT96ibgu2j9puU0aqItOCrpRNLnnzEvZ03WOJafeKKg0t6Eq5WMzhGNq+29YS2/zIZq6uf7VNGSlPoaNclHIhE2GsxXwvEA6tGrTSm31UqWkLXSkXuHbGtWw6uskSkzDRm31UmdKCrlQ5+nn/z3T8oKMltvfJvVwceHEheyh14bSgK1UORIQqL1p7NIdfN5x373zXEtObfVRZ0j50pcpY4zca5yvmEibOYq595aq86J2iSpWR5buX0+OjHpbY4acP06C6dQGK3P3m2oeuSkrvFFWqHGVLNl4vWtftHNNxDC/3eNmmjFRlpV0uSpWCb6RvvmIuYZKvmBc2IVfu77UrRpWWdrkodQG+3vY1/T7vZ4kljkmkdtXa591Xu1xUaWiXi1JlJDM7E59IH0tsYreJjOs8zqaMlDpHC7pSxZR3Ei24sLlXcg9V1GGLqixpl4tS5/HJxk/418J/WWIpY1Oo4VfDpoxUZaZdLkpdgDOZZ/Cf6G+JTes1jVE3jLIpI6WKVqKCbowZCgwWkS6FbG8D/A40F5E9pU1OKbuUVfeKUq5U7GGLxpj2wFtFbK8CvIe2+lUF9k70O/mKedrzaVrMVYVQrOJrjOkGLAS2FfGykcCVZZGUUq52Mv0k1SdXt8Rm953N4DaDbcpIqZIrbmu6MzAIaAt0ybvRGNMYmAg8B0wvq+SUcgXtXlGeorhdLpEi8nUR26cBS4BvS5+SUq7x6tpX8xXz9PHpWsxVhVWsFrqIZBe2zRjTBwgFrgKqFXUcY8xwYDjAxRfrfNDKHsmnk6n9svWOzgX3LKB/y/42ZaRU2SjVBUxjTHUcXSzPiMhRY0xwUa8XkfdwXDglJCREm0HK5bR7RXmy0k7O9RKwQ0RmlUUySpWX8T+Oz1fMMydkajFXHqW0Qwz7Ac2MMXn/V8QaY+aIyJBSHl+pUok/GU/91+pbYt/d/x09L+9pU0ZKlZ/SFvTegG+u5xfhuDB6O7C5lMdWqlTytsjrVq1LwpgEm7JRqvyVqqCLyNbcz40xyWe/3Soi+0pzbKUu1GNLHmP679bRs9kvZDvnH1fKU+ldncpjxKXE0eQ/TSyx1UNW07lZZ5syUsq1SlTQRST8PNv3ANoMUi6Xt3ulRd0WbB+13aZslLKHttBVhfavL//FJ5s+scS0e0VVVlrQVYW0K3EXl0+73BKLHhbN9Rddb1NGStlPC7qqcPJ2r3Ro0oF1D62zKRul3IcWdFVh3P7p7Sz5e4klpjcGKXVOae8UVarcbTm6BRNhrMX8LSAcjDEYYwgPD7cpO6Xch7bQlVvL271yR4s7WDRwEYQ5irlda+Iq5Y60oCu31OG/HfjlwC+WmHavKFU0LejKrUQfjKbd++0ssV2P7+LS2pfme21YWJir0lKqQtCCrtyCiFDlReslnfuvuZ+P7/q40H2031wpKy3oynYtprXg78S/LTHtXlGq5LSgK9us3rua0NmhltiBpw7QuGZjmzJSqmLTgq5crqDulVHtRjGt9zSbMlLKM2hBVy5V95W6JKYlWmLavaJU2dCCrlzi+53f0+uTXpbY0WeOUi+gnk0ZKeV5tKCrcpWVnYV3pPVjNr7zeCK7RdqUkVKeSwu6Kjd57/IE7V5RqjzpXC6qzC3YuiBfMU9+NlmLuVLlTFvoqsxkZGXg+5KvJfZqj1d5puMzNmWkVOWiBV2VCe1eUcp+2uWiSmV2zOx8xTz1uVQt5krZQFvo6oKkZaRRbVI1S+yd299hRMgImzJSSmlBVyWm3StKuSftclHFNu3XafmK+ennT2sxV8pNaAtdndeJMyeoOaWmJfbJXZ9w3zX32ZSRUqogWtBVkbR7RamKQ7tcVIEmrp6Yr5hnTMjQYq6UG9MWurJITEuk7it1LbGv7v2Kvlf2tSkjpVRxaUFXTnlb5H5efpwef9qmbJRSJaVdLopnf3g2XzHPeiFLi7lSFYy20Cuxw6mHafR6I0ts+QPL6X5pd5syUkqVhhb0Sipvi7xxjcYc+L8DNmWjlCoLWtArmeGLhvP+H+9bYtkvZGNM/uGJSqmKpUR96MaYocaYlXliDYwxC40xqcaYNGPMXGNMQJlmqUpt3/F9mAhjKebrHlyHhIkWc6U8RLFb6MaY9sBbwG95Ni0AGgPPAbWBCUA88HgZ5ahKKW/3yrUNrmXDvzfYlI1SqrwUq6AbY7oBC4FteeLdgVZASxE5dDbWEOiHFnTb3TP/HhZsXWCJ6Y1BSnmu4rbQOwODgLZAl1zxaKBDTjE/6xjgUybZqQuy49gOrph+hSUWMyKG1g1b25SRUsoVilvQI0Uk2xjTNndQRI4Dx/O8tgfwS1kkp0oub/dK1+Cu/Dj4R5uyUUq5UrEuiopIdnFeZ4y5FWgPTC1k+3BjTLQxJjo+Pr74Wap8wsPDLc9v+fCWfMVcwkSLuVKViBEpfp+qMSYc6CIiXQrYVhXYCPwtIr3Pd6yQkBCJjo4ufqbKwhiDiLDxyEZav2PtStk2chtXBF1RyJ5KqYrMGLNeREIK2laW49BfAeoCXcvwmKoIeVvkd111F18M+MKmbJRSdiuTuVyMMXcBo4BhIqK3G5aT8PBwjDGY+w2EW7dJmGgxV6qSK3VBN8a0AeYA00VEK0o5uv+x+x2FvMW52J4n9uhQRKUUUMouF2OMD/A5kATMNcbk7tfZKCLppTm+chARqryY53fvYpDftZArpc4pbR/6NZxrL67Ns+0SYE8pj1+phYeHc7TdUWZEz3DGvKt4kzEhg3AJty8xpZRbKtEol7Kko1yKtvnoZq6ZcY0llvRsErX8a9mTkFLKLbhqlIsqAwV1r8zpN4dBrQfZlJFSqqLQFYvcyKCFg6zF/CQQDoPbDM53I5FSSuWlLXQ3sP7gekLet/4FlTI2hZr+NbGrS0wpVfFoQbdRVnYW3pHWf4IF9yygf8v+NmWklKrItKDbpN9n/fh6+9fO55fXuZy/H/vb8pqwsDBXp6WUqsC0oLvY2n1r6TSrkyV2atwpqvpUzfda7TdXSpWEFnQXyczOxCfSOk38t/d9S+/m553HTCmlikULugt0ndOVlXtWOp+HXBTC78N+ty8hpZRH0oJejv63+3/c8tEtltiZ8Wfw9fK1KSOllCfTgl4OzmSewX+ivyX246Af6XqJziyslCo/WtDL2HXvXsefh/90Pr/l0lv44YEfbMxIKVVZaEEvI4t3LObOuXdaYhkTMvCuoqdYKeUaWm1K6VTGKQImBVhi6x5cR4emHWzKSClVWWlBL4XLpl7G7qTdzuf/uPIffHnvlzZmpJSqzLSgX4D5W+YzYMEASyzrhSyqGJ3rTCllHy3oJXDizAlqTqlpif0x/A/aNmprU0ZKKXWOFvRiqvNyHZJOJzmfD249mNn9ZtuXkFJK5aEF/TzmxMxhyNdDLLHsF7IxxtiTkFJKFUILeiGS0pKo80odS2zLo1toWa+lTRkppVTRtKAXoEpEFYRzC0uMajeKab2n2ZiRUkqdnxb0XN7+/W1GLhlpiWn3ilKqotCCDsSfjKf+a/UtsZ2P7eSyOpfZlJFSSpVcpS/oJsLa+n6u03NM6j7JpmyUUurCVdqC/sraV3h2+bOWmITpgsxKqYqr0hX0uJQ4mvyniSW278l9NA1salNGSilVNjziXvXirr1pIoylmE/sNhEJEy3mSimPYETs6WYICQmR6OjoMjmWMYaifo6wFWG8uPpFS0y7V5RSFZExZr2IhBS0zaO7XGKTYrl06qWW2KGnD9GwekObMlJKqfJTYbtcwsPDMcY4x4jnfJ/T/WIijKWYR90WhYSJFnOllMeq0AVdRJxdLTnfn+hwIt9QRAkTnrjxCTvSVEopl6nQBT237QnbMRGGN355wxlLGJ2gfeVKqUqjwl4UzbkQKiJUedH6e2nmnTN56LqHSpuiUkq5naIuipaohW6MGWqMWZkn5meMedsYk2iM2WGM6VWKXEvkm+3fWIp5gE8AEiZazJVSlVKxC7oxpj3wVgGbpgIDgceAKGCBMeaKMskuD+eF0OoGwqHvZ32d25KfTSZ1XGp5vK1SSlUIxepyMcZ0AxYCu4AUEelyNt4E2AsMEpFPzsZmAtkiMryoY15ol0vKmRQCpwQ6n2/890auaXBNiY+jlFIVUVl0uXQGBgHf5ImHAllA7qXuFwPdS5pkcfl7+9P/qv7wvWP0ihZzpZRyKO6NRZEikm2Mybsa8kXAbhFJyxXbDzQzxniJSFaZZJmLr5cvCwYsIHxreFkfWimlKrRitdBFJLuQTf5Acp5YGuAF1Mr7YmPMcGNMtDEmOj4+vgRp5lfc+VuUUqqyKO049DM4ulxySz/7tWreF4vIeyISIiIh9erVK+VbK6WUyq20Bf0ojm6X3Gqf/XqqlMdWSilVAqUt6BuBi40xuYt6W+A0kFTKYyullCqB0hb0GOAA8BSAMcYbGAb8KHbdgqqUUpVUqabPPTvyZTQw1xjTAmgEXIdjmKNSSikXKvXkXCIyD7gTx6iWU8CtIrKutMdVSilVMiVqoYtIeCHxJcCSskhIKaXUhbFttkVjzAlguy1vXjEFAQl2J1HB6DkrGT1fJWPX+WomIgWO+7ZzCbrthc1HoPIzxkTr+SoZPWclo+erZNzxfFXYBS6UUkpZaUFXSikPYWdBf8/G966I9HyVnJ6zktHzVTJud75suyiqlFKqbGmXi1JKeQgt6Eop5SFcVtDdbYFpd1fI+WpgjFlojEk1xqQZY+YaYwJsStGtFHS+8mxvY4zJMMYEuy4r91XU+TLG+BtjdhpjJrs4LbdVyP/HtsaYdcaYk8aY/caYSGOMrY1kl7y5OywwXZEUcb4WAK2B54DJwN1nv1ZqRZyvnO1VcFzAsvO+C7dxvvMFPAsEAJNck5F7K+h8GWOq4Vh68w+gF45z9X/AEy5PMJdy/4DnWmB6W554E+BhrAtMXwc8DRS5wLQnK+J8dQdaAS1F5NDZWEOgH/C4i9N0G4WdrzxGAle6JiP3dr7zZYy5BBgLjBKRE67MzR0Vcb46AHWAJ84utbnaGNMS6AP8x7VZnuOKFrrbLDBdQRR2vqKBDjnF/KxjgI+rEnNThZ0vAIwxjYGJOP6qUec5X8CbOIrXLJdl5N4KO19BOOqnyRXzx7GKm21cUdAjReTrAuJFLjDtgrzcVYHnS0SOi0jeVkIP4BfXpOW2Cvt85ZiGY+K4b12Uj7sr9HwZY3rimDk1AfjYGDNBr9EUer7W4lg7+QVjTA1jTFfgn8CnLs0uj3LvcinFAtPHyi8r91XE+bIwxtwKtAe6lW9G7q2o82WM6YPjL8GrgGouS8qNnefzlXM9pgGO2jAA6G+M6ZCn4VVpFHa+ROSAMeZh4BNgwtnwOyLyocuSK4CdV2RLtMC0OscYUxXHRZrvRGSF3fm4I2NMdWA68IyIHLU7H3dnjGkHtAHeF5FrRaQrjsZCK2Cwnbm5I2NMEPAqjoEK9wJhwGBjzDN25mXnVX9dYPrCvQLUBbranYgbewnYISLaF1w8zc9+fSMnICKrjTF/4xhZpayG4qhhA3KW2zTGxAOvGGPeFJEMO5Kys4WuC0xfAGPMXcAoYJiIHLA7HzfWD+hujBFjjACxZ+OxxpjZtmXlvlLPfo3NEz/Nub+c1TmXA9vyrJ28CaiO44KpLews6DHoAtMlYoxpA8wBpovIFzan4+5642gg5DxuPxu/HXjBrqTcWDQg5GqNG2PqA1cAv9uVlBuLB9rmGcDRE0dXsm0NUtu6XHSB6ZIxxvgAn+P4sMw1xuSeWH+jiGgrKhcR2Zr7uTEm+ey3W0Vkn+szcm8ictAYMxf4yBgzFkdhmgAcwdFPrKyWAOOAFcaYtTi6rPoBM0XktF1J2XrnnIjMM8ak4rgzTReYLto1QIuz36/Ns+0SYI9Ls1Ge6EEgEscd3HVwtNp72Vmg3JWIrDPG9MVxh+jTwEkcQxZH25mXTp+rlFIeQmdbVEopD6EFXSmlPIQWdKWU8hBa0JVSykNoQVdKKQ+hBV0ppTyEFnSllPIQ/w/gaYhCn0p73wAAAABJRU5ErkJggg==\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.988\n",
      "Model:                            OLS   Adj. R-squared:                  0.987\n",
      "Method:                 Least Squares   F-statistic:                     1291.\n",
      "Date:                Mon, 16 May 2022   Prob (F-statistic):           9.97e-17\n",
      "Time:                        23:03:49   Log-Likelihood:                -1.2941\n",
      "No. Observations:                  18   AIC:                             6.588\n",
      "Df Residuals:                      16   BIC:                             8.369\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept     -0.3878      0.426     -0.909      0.377      -1.292       0.516\n",
      "x              1.0092      0.028     35.927      0.000       0.950       1.069\n",
      "==============================================================================\n",
      "Omnibus:                        5.104   Durbin-Watson:                   2.532\n",
      "Prob(Omnibus):                  0.078   Jarque-Bera (JB):                2.648\n",
      "Skew:                           0.789   Prob(JB):                        0.266\n",
      "Kurtosis:                       4.019   Cond. No.                         100.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "假设检验结果为： 1290.750434422176\n",
      "落在拒绝域，所以X与Y的线性方程是显著的\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\goodboy\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\scipy\\stats\\_stats_py.py:1477: UserWarning: kurtosistest only valid for n>=20 ... continuing anyway, n=18\n",
      "  warnings.warn(\"kurtosistest only valid for n>=20 ... continuing \"\n"
     ]
    }
   ],
   "source": [
    "#4.7\n",
    "\n",
    "#参考csdn文章链接：https://blog.csdn.net/weixin_45807161/article/details/121459735\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "from statsmodels.formula.api import ols\n",
    "x = [17.1,10.5,13.8,15.7,11.9,10.4,15.0,16.0,17.8,\n",
    "     15.8,15.1,12.1,18.4,17.1,16.7,16.5,15.1,15.1]\n",
    "y = [16.7,10.4,13.5,15.7,11.6,10.2,14.5,15.8,17.6,\n",
    "     15.2,14.8,11.9,18.9,16.7,16.6,15.9,15.1,14.5]\n",
    "plt.plot(x, y, '+k', label = \"原始数据点\")\n",
    "p = np.polyfit(x, y, deg=1)  #拟合一次多项式\n",
    "print(\"拟合的多项式为:{}*x + {}\".format(p[0], p[1]))\n",
    "plt.rc('font', size = 16);\n",
    "plt.rc('font', family = 'SimHei')\n",
    "plt.plot(x, np.polyval(p,x), 'g-', label = \"拟合的直线\")\n",
    "plt.legend() #显示多个标签\n",
    "plt.show()\n",
    "#显著性检测 方法一,利用ols\n",
    "data = {'x':x, 'y':y}\n",
    "model = ols('y~x', data).fit()\n",
    "print(model.summary())\n",
    "#显著性检测 方法二，手推公式\n",
    "mean_x = np.mean(x)\n",
    "mean_y = np.mean(y)\n",
    "num1 = num2 = num3 = 0\n",
    "for i in range(0, 18):\n",
    "    num1 = num1 + x[i]*y[i]\n",
    "    num2 = num2 + y[i]**2\n",
    "    num3 = num3 + y[i]\n",
    "result1 = p[0]*(num1 - 18*mean_x*mean_y)\n",
    "result2 = num2 - p[1]*num3 - p[0]*num1\n",
    "F = result1 / (result2 / 16)\n",
    "print(\"假设检验结果为：\",F)\n",
    "if(F > 5.32):\n",
    "    print(\"落在拒绝域，所以X与Y的线性方程是显著的\")\n",
    "\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}